ToyMCSampler.cxx

Go to the documentation of this file.

// @(#)root/roostats:$Id$
// Author: Sven Kreiss    June 2010
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

/** \class RooStats::NuisanceParametersSampler
    \ingroup Roostats

Helper class for ToyMCSampler. Handles all of the nuisance parameter related
functions. Once instantiated, it gives a new nuisance parameter point
at each call to nextPoint(...).
*/

/** \class RooStats::ToyMCSampler
    \ingroup Roostats

ToyMCSampler is an implementation of the TestStatSampler interface.
It generates Toy Monte Carlo for a given parameter point and evaluates a
TestStatistic.

For parallel runs, ToyMCSampler can be given an instance of ProofConfig
and then run in parallel using proof or proof-lite. Internally, it uses
ToyMCStudy with the RooStudyManager.
*/

#include "RooStats/ToyMCSampler.h"

#include "RooMsgService.h"

#include "RooDataHist.h"

#include "RooRealVar.h"

#include "TCanvas.h"
#include "RooPlot.h"
#include "RooRandom.h"

#include "RooStudyManager.h"
#include "RooStats/ToyMCStudy.h"
#include "RooStats/DetailedOutputAggregator.h"
#include "RooStats/RooStatsUtils.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"

#include "TMath.h"


using namespace RooFit;
using namespace std;


ClassImp(RooStats::ToyMCSampler);

namespace RooStats {

////////////////////////////////////////////////////////////////////////////////
/// Assigns new nuisance parameter point to members of nuisPoint.
/// nuisPoint can be more objects than just the nuisance
/// parameters.

void NuisanceParametersSampler::NextPoint(RooArgSet& nuisPoint, Double_t& weight) {

   // check whether to get new set of nuisanceParPoints
   if (fIndex >= fNToys) {
      Refresh();
      fIndex = 0;
   }

   // get value
   nuisPoint =  *fPoints->get(fIndex++);
   weight = fPoints->weight();

   // check whether result will have any influence
   if(fPoints->weight() == 0.0) {
      oocoutI((TObject*)NULL,Generation) << "Weight 0 encountered. Skipping." << endl;
      NextPoint(nuisPoint, weight);
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Creates the initial set of nuisance parameter points. It also refills the
/// set with new parameter points if called repeatedly. This helps with
/// adaptive sampling as the required number of nuisance parameter points
/// might increase during the run.

void NuisanceParametersSampler::Refresh() {

   if (!fPrior || !fParams) return;

   if (fPoints) delete fPoints;

   if (fExpected) {
      // UNDER CONSTRUCTION
      oocoutI((TObject*)NULL,InputArguments) << "Using expected nuisance parameters." << endl;

      int nBins = fNToys;

      // From FeldmanCousins.cxx:
      // set nbins for the POI
      TIter it2 = fParams->createIterator();
      RooRealVar *myarg2;
      while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
        myarg2->setBins(nBins);
      }


      fPoints = fPrior->generate(
         *fParams,
         AllBinned(),
         ExpectedData(),
         NumEvents(1) // for Asimov set, this is only a scale factor
      );
      if(fPoints->numEntries() != fNToys) {
         fNToys = fPoints->numEntries();
         oocoutI((TObject*)NULL,InputArguments) <<
            "Adjusted number of toys to number of bins of nuisance parameters: " << fNToys << endl;
      }

/*
      // check
      TCanvas *c1 = new TCanvas;
      RooPlot *p = dynamic_cast<RooRealVar*>(fParams->first())->frame();
      fPoints->plotOn(p);
      p->Draw();
      for(int x=0; x < fPoints->numEntries(); x++) {
         fPoints->get(x)->Print("v");
         cout << fPoints->weight() << endl;
      }
*/

   }else{
      oocoutI((TObject*)NULL,InputArguments) << "Using randomized nuisance parameters." << endl;

      fPoints = fPrior->generate(*fParams, fNToys);
   }
}

Bool_t ToyMCSampler::fgAlwaysUseMultiGen = kFALSE ;

////////////////////////////////////////////////////////////////////////////////

void ToyMCSampler::SetAlwaysUseMultiGen(Bool_t flag) { fgAlwaysUseMultiGen = flag ; }

////////////////////////////////////////////////////////////////////////////////
/// Proof constructor. Do not use.

ToyMCSampler::ToyMCSampler() : fSamplingDistName("SD"), fNToys(1)
{

   fPdf = NULL;
   fParametersForTestStat = NULL;
   fPriorNuisance = NULL;
   fNuisancePars = NULL;
   fObservables = NULL;
   fGlobalObservables = NULL;

   fSize = 0.05;
   fNEvents = 0;
   fGenerateBinned = kFALSE;
   fGenerateBinnedTag = "";
   fGenerateAutoBinned = kTRUE;
   fExpectedNuisancePar = kFALSE;

   fToysInTails = 0.0;
   fMaxToys = RooNumber::infinity();
   fAdaptiveLowLimit = -RooNumber::infinity();
   fAdaptiveHighLimit = RooNumber::infinity();

   fProtoData = NULL;

   fProofConfig = NULL;
   fNuisanceParametersSampler = NULL;

   _allVars = NULL ;
   _gs1 = NULL ;
   _gs2 = NULL ;
   _gs3 = NULL ;
   _gs4 = NULL ;

   //suppress messages for num integration of Roofit
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);

   fUseMultiGen = kFALSE ;
}

////////////////////////////////////////////////////////////////////////////////

ToyMCSampler::ToyMCSampler(TestStatistic &ts, Int_t ntoys) : fSamplingDistName(ts.GetVarName().Data()), fNToys(ntoys)
{
   fPdf = NULL;
   fParametersForTestStat = NULL;
   fPriorNuisance = NULL;
   fNuisancePars = NULL;
   fObservables = NULL;
   fGlobalObservables = NULL;

   fSize = 0.05;
   fNEvents = 0;
   fGenerateBinned = kFALSE;
   fGenerateBinnedTag = "";
   fGenerateAutoBinned = kTRUE;
   fExpectedNuisancePar = kFALSE;

   fToysInTails = 0.0;
   fMaxToys = RooNumber::infinity();
   fAdaptiveLowLimit = -RooNumber::infinity();
   fAdaptiveHighLimit = RooNumber::infinity();

   fProtoData = NULL;

   fProofConfig = NULL;
   fNuisanceParametersSampler = NULL;

   _allVars = NULL ;
   _gs1 = NULL ;
   _gs2 = NULL ;
   _gs3 = NULL ;
   _gs4 = NULL ;

   //suppress messages for num integration of Roofit
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);

   fUseMultiGen = kFALSE ;

   AddTestStatistic(&ts);
}

////////////////////////////////////////////////////////////////////////////////

ToyMCSampler::~ToyMCSampler() {
   if(fNuisanceParametersSampler) delete fNuisanceParametersSampler;

   ClearCache();
}

////////////////////////////////////////////////////////////////////////////////
/// only checks, no guessing/determination (do this in calculators,
/// e.g. using ModelConfig::GuessObsAndNuisance(...))

Bool_t ToyMCSampler::CheckConfig(void) {
   bool goodConfig = true;

   if(fTestStatistics.size() == 0 || fTestStatistics[0] == NULL) { ooccoutE((TObject*)NULL,InputArguments) << "Test statistic not set." << endl; goodConfig = false; }
   if(!fObservables) { ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl; goodConfig = false; }
   if(!fParametersForTestStat) { ooccoutE((TObject*)NULL,InputArguments) << "Parameter values used to evaluate the test statistic are not set." << endl; goodConfig = false; }
   if(!fPdf) { ooccoutE((TObject*)NULL,InputArguments) << "Pdf not set." << endl; goodConfig = false; }


   //ooccoutI((TObject*)NULL,InputArguments) << "ToyMCSampler configuration:" << endl;
   //ooccoutI((TObject*)NULL,InputArguments) << "Pdf from SetPdf: " << fPdf << endl;
   // for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
   //   ooccoutI((TObject*)NULL,InputArguments) << "test statistics["<<i<<"]: " << fTestStatistics[i] << endl;
   // }
   //ooccoutI((TObject*)NULL,InputArguments) << endl;

   return goodConfig;
}

////////////////////////////////////////////////////////////////////////////////
/// Evaluate all test statistics, returning result and any detailed output.
/// PDF parameter values are saved in case they are modified by
/// TestStatistic::Evaluate (eg. SimpleLikelihoodRatioTestStat).

RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi) {
   DetailedOutputAggregator detOutAgg;
   const RooArgList* allTS = EvaluateAllTestStatistics(data, poi, detOutAgg);
   if (!allTS) return 0;
   // no need to delete allTS, it is deleted in destructor of detOutAgg
   return  dynamic_cast<RooArgList*>(allTS->snapshot());
}

////////////////////////////////////////////////////////////////////////////////

const RooArgList* ToyMCSampler::EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg) {
   RooArgSet *allVars = fPdf ? fPdf->getVariables() : nullptr;
   RooArgSet *saveAll = allVars ? allVars->snapshot() : nullptr;
   for( unsigned int i = 0; i < fTestStatistics.size(); i++ ) {
      if( fTestStatistics[i] == NULL ) continue;
      TString name( TString::Format("%s_TS%u", fSamplingDistName.c_str(), i) );
      std::unique_ptr<RooArgSet> parForTS(poi.snapshot());
      RooRealVar ts( name, fTestStatistics[i]->GetVarName(), fTestStatistics[i]->Evaluate( data, *parForTS ) );
      RooArgList tset(ts);
      detOutAgg.AppendArgSet(&tset);
      if (const RooArgSet* detOut = fTestStatistics[i]->GetDetailedOutput()) {
        name.Append("_");
        detOutAgg.AppendArgSet(detOut, name);
      }
      if (saveAll) *allVars = *saveAll;  // restore values, perhaps modified by fTestStatistics[i]->Evaluate()
   }
   delete saveAll;
   delete allVars;
   return detOutAgg.GetAsArgList();
}

////////////////////////////////////////////////////////////////////////////////

SamplingDistribution* ToyMCSampler::GetSamplingDistribution(RooArgSet& paramPointIn) {
   if(fTestStatistics.size() > 1) {
      oocoutW((TObject*)NULL, InputArguments) << "Multiple test statistics defined, but only one distribution will be returned." << endl;
      for( unsigned int i=0; i < fTestStatistics.size(); i++ ) {
         oocoutW((TObject*)NULL, InputArguments) << " \t test statistic: " << fTestStatistics[i] << endl;
      }
   }

   RooDataSet* r = GetSamplingDistributions(paramPointIn);
   if(r == NULL || r->numEntries() == 0) {
      oocoutW((TObject*)NULL, Generation) << "no sampling distribution generated" << endl;
      return NULL;
   }

   SamplingDistribution* samp = new SamplingDistribution( r->GetName(), r->GetTitle(), *r );
   delete r;
   return samp;
}

////////////////////////////////////////////////////////////////////////////////
/// Use for serial and parallel runs.

RooDataSet* ToyMCSampler::GetSamplingDistributions(RooArgSet& paramPointIn)
{

   // ======= S I N G L E   R U N ? =======
   if(!fProofConfig)
      return GetSamplingDistributionsSingleWorker(paramPointIn);

   // ======= P A R A L L E L   R U N =======
   if (!CheckConfig()){
      oocoutE((TObject*)NULL, InputArguments)
         << "Bad COnfiguration in ToyMCSampler "
         << endl;
      return nullptr;
   }

   // turn adaptive sampling off if given
   if(fToysInTails) {
      fToysInTails = 0;
      oocoutW((TObject*)NULL, InputArguments)
         << "Adaptive sampling in ToyMCSampler is not supported for parallel runs."
         << endl;
   }

   // adjust number of toys on the slaves to keep the total number of toys constant
   Int_t totToys = fNToys;
   fNToys = (int)ceil((double)fNToys / (double)fProofConfig->GetNExperiments()); // round up

   // create the study instance for parallel processing
   ToyMCStudy* toymcstudy = new ToyMCStudy ;
   toymcstudy->SetToyMCSampler(*this);
   toymcstudy->SetParamPoint(paramPointIn);
   toymcstudy->SetRandomSeed(RooRandom::randomGenerator()->Integer(TMath::Limits<unsigned int>::Max() ) );

   // temporary workspace for proof to avoid messing with TRef
   RooWorkspace w(fProofConfig->GetWorkspace());
   RooStudyManager studymanager(w, *toymcstudy);
   studymanager.runProof(fProofConfig->GetNExperiments(), fProofConfig->GetHost(), fProofConfig->GetShowGui());

   RooDataSet* output = toymcstudy->merge();

   // reset the number of toys
   fNToys = totToys;

   delete toymcstudy;
   return output;
}

////////////////////////////////////////////////////////////////////////////////
/// This is the main function for serial runs. It is called automatically
/// from inside GetSamplingDistribution when no ProofConfig is given.
/// You should not call this function yourself. This function should
/// be used by ToyMCStudy on the workers (ie. when you explicitly want
/// a serial run although ProofConfig is present).
///

RooDataSet* ToyMCSampler::GetSamplingDistributionsSingleWorker(RooArgSet& paramPointIn)
{
  // Make sure the cache is clear. It is important to clear it here, because
  // the cache might be invalid even when just the firstPOI was changed, for which
  // no accessor has to be called. (Fixes a bug when ToyMCSampler is
  // used with the Neyman Construction)
   ClearCache();

   if (!CheckConfig()){
      oocoutE((TObject*)NULL, InputArguments)
         << "Bad COnfiguration in ToyMCSampler "
         << endl;
      return nullptr;
   }

   // important to cache the paramPoint b/c test statistic might
   // modify it from event to event
   RooArgSet *paramPoint = (RooArgSet*) paramPointIn.snapshot();
   RooArgSet *allVars = fPdf->getVariables();
   RooArgSet *saveAll = (RooArgSet*) allVars->snapshot();


   DetailedOutputAggregator detOutAgg;

   // counts the number of toys in the limits set for adaptive sampling
   // (taking weights into account; always on first test statistic)
   Double_t toysInTails = 0.0;

   for (Int_t i = 0; i < fMaxToys; ++i) {
      // need to check at the beginning for case that zero toys are requested
      if (toysInTails >= fToysInTails  &&  i+1 > fNToys) break;

      // status update
      if ( i% 500 == 0 && i>0 ) {
         oocoutP((TObject*)0,Generation) << "generated toys: " << i << " / " << fNToys;
         if (fToysInTails) ooccoutP((TObject*)0,Generation) << " (tails: " << toysInTails << " / " << fToysInTails << ")" << std::endl;
         else ooccoutP((TObject*)0,Generation) << endl;
      }

      // TODO: change this treatment to keep track of all values so that the threshold
      // for adaptive sampling is counted for all distributions and not just the
      // first one.
      Double_t valueFirst = -999.0, weight = 1.0;

      // set variables to requested parameter point
      *allVars = *saveAll; // important for example for SimpleLikelihoodRatioTestStat

      RooAbsData* toydata = GenerateToyData(*paramPoint, weight);
      if (i == 0 && !fPdf->canBeExtended() && dynamic_cast<RooSimultaneous*>(fPdf)) {
        const RooArgSet* toySet = toydata->get();
        if (std::none_of(toySet->begin(), toySet->end(), [](const RooAbsArg* arg){
          return dynamic_cast<const RooAbsCategory*>(arg) != nullptr;
        }))
          oocoutE((TObject*)nullptr, Generation) << "ToyMCSampler: Generated toy data didn't contain a category variable, although"
            " a simultaneous PDF is in use. To generate events for a simultaneous PDF, all components need to be"
            " extended. Otherwise, the number of events to generate per component cannot be determined." << std::endl;
      }

      *allVars = *fParametersForTestStat;

      const RooArgList* allTS = EvaluateAllTestStatistics(*toydata, *fParametersForTestStat, detOutAgg);
      if (allTS->getSize() > Int_t(fTestStatistics.size()))
        detOutAgg.AppendArgSet( fGlobalObservables, "globObs_" );
      if (RooRealVar* firstTS = dynamic_cast<RooRealVar*>(allTS->first()))
         valueFirst = firstTS->getVal();

      delete toydata;

      // check for nan
      if(valueFirst != valueFirst) {
         oocoutW((TObject*)NULL, Generation) << "skip: " << valueFirst << ", " << weight << endl;
         continue;
      }

      detOutAgg.CommitSet(weight);

      // adaptive sampling checks
      if (valueFirst <= fAdaptiveLowLimit  ||  valueFirst >= fAdaptiveHighLimit) {
         if(weight >= 0.) toysInTails += weight;
         else toysInTails += 1.;
      }
   }

   // clean up
   *allVars = *saveAll;
   delete saveAll;
   delete allVars;
   delete paramPoint;

   return detOutAgg.GetAsDataSet(fSamplingDistName, fSamplingDistName);
}

////////////////////////////////////////////////////////////////////////////////

void ToyMCSampler::GenerateGlobalObservables(RooAbsPdf& pdf) const {


   if(!fGlobalObservables  ||  fGlobalObservables->getSize()==0) {
      ooccoutE((TObject*)NULL,InputArguments) << "Global Observables not set." << endl;
      return;
   }


   if (fUseMultiGen || fgAlwaysUseMultiGen) {

      // generate one set of global observables and assign it
      // has problem for sim pdfs
      RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>( &pdf );
      if (!simPdf) {
         RooDataSet *one = pdf.generate(*fGlobalObservables, 1);

         const RooArgSet *values = one->get(0);
         if (!_allVars) {
            _allVars = pdf.getVariables();
         }
         *_allVars = *values;
         delete one;

      } else {

         if (_pdfList.size() == 0) {
            RooCategory& channelCat = (RooCategory&)simPdf->indexCat();
            int nCat = channelCat.numTypes();
            for (int i=0; i < nCat; ++i){
               channelCat.setIndex(i);
               RooAbsPdf* pdftmp = simPdf->getPdf(channelCat.getLabel());
               assert(pdftmp);
               RooArgSet* globtmp = pdftmp->getObservables(*fGlobalObservables);
               RooAbsPdf::GenSpec* gs = pdftmp->prepareMultiGen(*globtmp, NumEvents(1));
               _pdfList.push_back(pdftmp);
               _obsList.push_back(globtmp);
               _gsList.push_back(gs);
            }
         }

         list<RooArgSet*>::iterator oiter = _obsList.begin();
         list<RooAbsPdf::GenSpec*>::iterator giter = _gsList.begin();
         for (list<RooAbsPdf*>::iterator iter = _pdfList.begin(); iter != _pdfList.end(); ++iter, ++giter, ++oiter) {
            //RooDataSet* tmp = (*iter)->generate(**oiter,1) ;
            RooDataSet* tmp = (*iter)->generate(**giter);
            **oiter = *tmp->get(0);
            delete tmp;
         }
      }


   } else {

      // not using multigen for global observables
      RooDataSet* one = pdf.generateSimGlobal( *fGlobalObservables, 1 );
      const RooArgSet *values = one->get(0);
      RooArgSet* allVars = pdf.getVariables();
      *allVars = *values;
      delete allVars;
      delete one;

   }
}

////////////////////////////////////////////////////////////////////////////////
/// This method generates a toy data set for the given parameter point taking
/// global observables into account.
/// The values of the generated global observables remain in the pdf's variables.
/// They have to have those values for the subsequent evaluation of the
/// test statistics.

RooAbsData* ToyMCSampler::GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const {

   if(!fObservables) {
      ooccoutE((TObject*)NULL,InputArguments) << "Observables not set." << endl;
      return NULL;
   }

   // assign input paramPoint
   RooArgSet* allVars = fPdf->getVariables();
   *allVars = paramPoint;


   // create nuisance parameter points
   if(!fNuisanceParametersSampler && fPriorNuisance && fNuisancePars) {
      fNuisanceParametersSampler = new NuisanceParametersSampler(fPriorNuisance, fNuisancePars, fNToys, fExpectedNuisancePar);
      if ((fUseMultiGen || fgAlwaysUseMultiGen) &&  fNuisanceParametersSampler )
         oocoutI((TObject*)NULL,InputArguments) << "Cannot use multigen when nuisance parameters vary for every toy" << endl;
   }

   // generate global observables
   RooArgSet observables(*fObservables);
   if(fGlobalObservables  &&  fGlobalObservables->getSize()) {
      observables.remove(*fGlobalObservables);
      GenerateGlobalObservables(pdf);
   }

   // save values to restore later.
   // but this must remain after(!) generating global observables
   const RooArgSet* saveVars = (const RooArgSet*)allVars->snapshot();

   if(fNuisanceParametersSampler) { // use nuisance parameters?
      // Construct a set of nuisance parameters that has the parameters
      // in the input paramPoint removed. Therefore, no parameter in
      // paramPoint is randomized.
      // Therefore when a parameter is given (should be held fixed),
      // but is also in the list of nuisance parameters, the parameter
      // will be held fixed. This is useful for debugging to hold single
      // parameters fixed although under "normal" circumstances it is
      // randomized.
      RooArgSet allVarsMinusParamPoint(*allVars);
      allVarsMinusParamPoint.remove(paramPoint, kFALSE, kTRUE); // match by name

      // get nuisance parameter point and weight
      fNuisanceParametersSampler->NextPoint(allVarsMinusParamPoint, weight);


   }else{
      weight = 1.0;
   }

   RooAbsData *data = Generate(pdf, observables);

   // We generated the data with the randomized nuisance parameter (if hybrid)
   // but now we are setting the nuisance parameters back to where they were.
   *allVars = *saveVars;
   delete allVars;
   delete saveVars;

   return data;
}

////////////////////////////////////////////////////////////////////////////////
/// This is the generate function to use in the context of the ToyMCSampler
/// instead of the standard RooAbsPdf::generate(...).
/// It takes into account whether the number of events is given explicitly
/// or whether it should use the expected number of events. It also takes
/// into account the option to generate a binned data set (*i.e.* RooDataHist).

RooAbsData* ToyMCSampler::Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet* protoData, int forceEvents) const {

  if(fProtoData) {
    protoData = fProtoData;
    forceEvents = protoData->numEntries();
  }

  RooAbsData *data = NULL;
  int events = forceEvents;
  if(events == 0) events = fNEvents;

  // cannot use multigen when the nuisance parameters change for every toy
  bool useMultiGen = (fUseMultiGen || fgAlwaysUseMultiGen) && !fNuisanceParametersSampler;

  if (events == 0) {
    if (pdf.canBeExtended() && pdf.expectedEvents(observables) > 0) {
      if(fGenerateBinned) {
        if(protoData) data = pdf.generate(observables, AllBinned(), Extended(), ProtoData(*protoData, true, true));
        else          data = pdf.generate(observables, AllBinned(), Extended());
      } else {
        if (protoData) {
          if (useMultiGen) {
            if (!_gs2) { _gs2 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)) ; }
            data = pdf.generate(*_gs2) ;
          } else {
            data = pdf.generate                    (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
          }
        } else {
          if (useMultiGen) {
            if (!_gs1) { _gs1 = pdf.prepareMultiGen(observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) ) ; }
            data = pdf.generate(*_gs1) ;
          } else {
            data = pdf.generate                    (observables, Extended(), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag) );
          }

        }
      }
    } else {
      oocoutE((TObject*)0,InputArguments)
                << "ToyMCSampler: Error : pdf is not extended and number of events per toy is zero"
                << endl;
    }
  } else {
    if (fGenerateBinned) {
      if(protoData) data = pdf.generate(observables, events, AllBinned(), ProtoData(*protoData, true, true));
      else          data = pdf.generate(observables, events, AllBinned());
    } else {
      if (protoData) {
        if (useMultiGen) {
          if (!_gs3) { _gs3 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true)); }
          data = pdf.generate(*_gs3) ;
        } else {
          data = pdf.generate                    (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag), ProtoData(*protoData, true, true));
        }
      } else {
        if (useMultiGen) {
          if (!_gs4) { _gs4 = pdf.prepareMultiGen(observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag)); }
          data = pdf.generate(*_gs4) ;
        } else {
          data = pdf.generate                    (observables, NumEvents(events), AutoBinned(fGenerateAutoBinned), GenBinned(fGenerateBinnedTag));
        }
      }
    }
  }

  // in case of number counting print observables
  // if (data->numEntries() == 1) {
  //    std::cout << "generate observables : ";
  //    RooStats::PrintListContent(*data->get(0), std::cout);
  // }

  return data;
}

////////////////////////////////////////////////////////////////////////////////
/// Extended interface to append to sampling distribution more samples

SamplingDistribution* ToyMCSampler::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;
}

////////////////////////////////////////////////////////////////////////////////
/// clear the cache obtained from the pdf used for speeding the toy and global observables generation
/// needs to be called every time the model pdf (fPdf) changes

void ToyMCSampler::ClearCache() {

   if (_gs1) delete _gs1;
   _gs1 = 0;
   if (_gs2) delete _gs2;
   _gs2 = 0;
   if (_gs3) delete _gs3;
   _gs3 = 0;
   if (_gs4) delete _gs4;
   _gs4 = 0;

   // no need to delete the _pdfList since it is managed by the RooSimultaneous object
   if (_pdfList.size() > 0) {
      std::list<RooArgSet*>::iterator oiter = _obsList.begin();
      for (std::list<RooAbsPdf::GenSpec*>::iterator giter  = _gsList.begin(); giter != _gsList.end(); ++giter, ++oiter) {
         delete *oiter;
         delete *giter;
      }
      _pdfList.clear();
      _obsList.clear();
      _gsList.clear();
   }

   //LM: is this set really needed ??
   if (_allVars) delete _allVars;
   _allVars = 0;

}

} // end namespace RooStats