// @(#)root/roostats:$Id$
// Author: Sven Kreiss and Kyle Cranmer    June 2010
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
// Additions and modifications by Mario Pelliccioni
/*************************************************************************
 * 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.             *
 *************************************************************************/

#ifndef ROOSTATS_ToyMCSampler
#define ROOSTATS_ToyMCSampler

//_________________________________________________
/*
BEGIN_HTML
<p>
ToyMCSampler is an implementation of the TestStatSampler interface.
It generates Toy Monte Carlo for a given parameter point and evaluates a
TestStatistic.
</p>

<p>
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.
</p>
END_HTML
*/
//

#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif

#include <vector>
#include <sstream>

#include "RooStats/TestStatSampler.h"
#include "RooStats/SamplingDistribution.h"
#include "RooStats/TestStatistic.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/ProofConfig.h"

#include "RooWorkspace.h"
#include "RooMsgService.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"

#include "RooDataSet.h"

namespace RooStats {

  class DetailedOutputAggregator;

// only used inside ToyMCSampler, ie "private" in the cxx file
class NuisanceParametersSampler {
   // Helper for ToyMCSampler. Handles all of the nuisance parameter related
   // functions. Once instantiated, it gives a new nuisance parameter point
   // at each call to nextPoint(...).

   public:
      NuisanceParametersSampler(RooAbsPdf *prior=NULL, const RooArgSet *parameters=NULL, Int_t nToys=1000, Bool_t asimov=kFALSE) :
         fPrior(prior),
         fParams(parameters),
         fNToys(nToys),
         fExpected(asimov),
         fPoints(NULL),
         fIndex(0)
      {
         if(prior) Refresh();
      }
      virtual ~NuisanceParametersSampler() {
         if(fPoints) { delete fPoints; fPoints = NULL; }
      }

      void NextPoint(RooArgSet& nuisPoint, Double_t& weight);

   protected:
      void Refresh();

   private:
      RooAbsPdf *fPrior;           // prior for nuisance parameters
      const RooArgSet *fParams;    // nuisance parameters
      Int_t fNToys;
      Bool_t fExpected;

      RooAbsData *fPoints;         // generated nuisance parameter points
      Int_t fIndex;                // current index in fPoints array
};

class ToyMCSampler: public TestStatSampler {

   public:

      ToyMCSampler();
      ToyMCSampler(TestStatistic &ts, Int_t ntoys);
      virtual ~ToyMCSampler();
   
      static void SetAlwaysUseMultiGen(Bool_t flag);

      void SetUseMultiGen(Bool_t flag) { fUseMultiGen = flag ; }

      // main interface
      virtual SamplingDistribution* GetSamplingDistribution(RooArgSet& paramPoint);
      virtual RooDataSet* GetSamplingDistributions(RooArgSet& paramPoint);
      virtual RooDataSet* GetSamplingDistributionsSingleWorker(RooArgSet& paramPoint);

      virtual SamplingDistribution* AppendSamplingDistribution(
         RooArgSet& allParameters, 
         SamplingDistribution* last, 
			Int_t additionalMC
		);


      // The pdf can be NULL in which case the density from SetPdf()
      // is used. The snapshot and TestStatistic is also optional.
      virtual void AddTestStatistic(TestStatistic* t = NULL) {
         if( t == NULL ) {
            oocoutI((TObject*)0,InputArguments) << "No test statistic given. Doing nothing." << std::endl;
            return;
         }

         //if( t == NULL && fTestStatistics.size() >= 1 ) t = fTestStatistics[0];
         
         fTestStatistics.push_back( t );
      }      



      // generates toy data
      //   without weight
      virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, RooAbsPdf& pdf) const {
         if(fExpectedNuisancePar) oocoutE((TObject*)NULL,InputArguments) << "ToyMCSampler: using expected nuisance parameters but ignoring weight. Use GetSamplingDistribution(paramPoint, weight) instead." << std::endl;
         double weight;
         return GenerateToyData(paramPoint, weight, pdf);
      }
      virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint) const { return GenerateToyData(paramPoint,*fPdf); }
      //   with weight
      virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const;
      virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const { return GenerateToyData(paramPoint,weight,*fPdf); }

      // generate global observables
      virtual void GenerateGlobalObservables(RooAbsPdf& pdf) const;


      // Main interface to evaluate the test statistic on a dataset
      virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI, int i ) {
         return fTestStatistics[i]->Evaluate(data, nullPOI);
      }
      virtual Double_t EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI) { return EvaluateTestStatistic( data,nullPOI, 0 ); }
      virtual RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi);


      virtual TestStatistic* GetTestStatistic(unsigned int i) const {
         if( fTestStatistics.size() <= i ) return NULL;
         return fTestStatistics[i];
      }
      virtual TestStatistic* GetTestStatistic(void) const { return GetTestStatistic(0); }
      
      virtual Double_t ConfidenceLevel() const { return 1. - fSize; }
      virtual void Initialize(
         RooAbsArg& /*testStatistic*/,
         RooArgSet& /*paramsOfInterest*/,
         RooArgSet& /*nuisanceParameters*/
      ) {}

      virtual Int_t GetNToys(void) { return fNToys; }
      virtual void SetNToys(const Int_t ntoy) { fNToys = ntoy; }
      virtual void SetNEventsPerToy(const Int_t nevents) {
         // Forces n events even for extended PDFs. Set NEvents=0 to
         // use the Poisson distributed events from the extended PDF.
         fNEvents = nevents;
      }


      // Set the Pdf, add to the the workspace if not already there
      virtual void SetParametersForTestStat(const RooArgSet& nullpoi) {
         if( fParametersForTestStat ) delete fParametersForTestStat;
         fParametersForTestStat = (const RooArgSet*)nullpoi.snapshot();
      }

      virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; ClearCache(); }

      // How to randomize the prior. Set to NULL to deactivate randomization.
      virtual void SetPriorNuisance(RooAbsPdf* pdf) { 
         fPriorNuisance = pdf; 
         if (fNuisanceParametersSampler) { 
            delete fNuisanceParametersSampler; 
            fNuisanceParametersSampler = NULL; 
         }
      }
      // specify the nuisance parameters (eg. the rest of the parameters)
      virtual void SetNuisanceParameters(const RooArgSet& np) { fNuisancePars = &np; }
      // specify the observables in the dataset (needed to evaluate the test statistic)
      virtual void SetObservables(const RooArgSet& o) { fObservables = &o; }
      // specify the conditional observables
      virtual void SetGlobalObservables(const RooArgSet& o) { fGlobalObservables = &o; }


      // set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
      virtual void SetTestSize(Double_t size) { fSize = size; }
      // set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
      virtual void SetConfidenceLevel(Double_t cl) { fSize = 1. - cl; }

      // Set the TestStatistic (want the argument to be a function of the data & parameter points
      virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i) {
         if( fTestStatistics.size() < i ) {
            oocoutE((TObject*)NULL,InputArguments) << "Cannot set test statistic for this index." << std::endl;
            return;
         }
	 if( fTestStatistics.size() == i)
		 fTestStatistics.push_back(testStatistic);
	 else
		 fTestStatistics[i] = testStatistic;
      }
      virtual void SetTestStatistic(TestStatistic *t) { return SetTestStatistic(t,0); }

      virtual void SetExpectedNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }
      virtual void SetAsimovNuisancePar(Bool_t i = kTRUE) { fExpectedNuisancePar = i; }

      // Checks for sufficient information to do a GetSamplingDistribution(...).
      Bool_t CheckConfig(void);

      // control to use bin data generation (=> see RooFit::AllBinned() option)
      void SetGenerateBinned(bool binned = true) { fGenerateBinned = binned; }
      // name of the tag for individual components to be generated binned (=> see RooFit::GenBinned() option)
      void SetGenerateBinnedTag( const char* binnedTag = "" ) { fGenerateBinnedTag = binnedTag; }
      // set auto binned generation (=> see RooFit::AutoBinned() option)
      void SetGenerateAutoBinned( Bool_t autoBinned = kTRUE ) { fGenerateAutoBinned = autoBinned; }

      // Set the name of the sampling distribution used for plotting
      void SetSamplingDistName(const char* name) { if(name) fSamplingDistName = name; }
      std::string GetSamplingDistName(void) { return fSamplingDistName; }

      // This option forces a maximum number of total toys.
      void SetMaxToys(Double_t t) { fMaxToys = t; }

      void SetToysLeftTail(Double_t toys, Double_t threshold) {
         fToysInTails = toys;
         fAdaptiveLowLimit = threshold;
         fAdaptiveHighLimit = RooNumber::infinity();
      }
      void SetToysRightTail(Double_t toys, Double_t threshold) {
         fToysInTails = toys;
         fAdaptiveHighLimit = threshold;
         fAdaptiveLowLimit = -RooNumber::infinity();
      }
      void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold) {
         fToysInTails = toys;
         fAdaptiveHighLimit = high_threshold;
         fAdaptiveLowLimit = low_threshold;
      }

      // calling with argument or NULL deactivates proof
      void SetProofConfig(ProofConfig *pc = NULL) { fProofConfig = pc; }

      void SetProtoData(const RooDataSet* d) { fProtoData = d; }
      
   protected:

      const RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg);

      // helper for GenerateToyData
      RooAbsData* Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=NULL, int forceEvents=0) const;

      // helper method for clearing  the cache
      virtual void ClearCache();


      // densities, snapshots, and test statistics to reweight to
      RooAbsPdf *fPdf; // model (can be alt or null)
      const RooArgSet* fParametersForTestStat;
      std::vector<TestStatistic*> fTestStatistics;

      std::string fSamplingDistName; // name of the model
      RooAbsPdf *fPriorNuisance; // prior pdf for nuisance parameters
      const RooArgSet *fNuisancePars;
      const RooArgSet *fObservables;
      const RooArgSet *fGlobalObservables;
      Int_t fNToys; // number of toys to generate
      Int_t fNEvents; // number of events per toy (may be ignored depending on settings)
      Double_t fSize;
      Bool_t fExpectedNuisancePar; // whether to use expectation values for nuisance parameters (ie Asimov data set)
      Bool_t fGenerateBinned;
      TString fGenerateBinnedTag;
      Bool_t fGenerateAutoBinned;

      // minimum no of toys in tails for adaptive sampling
      // (taking weights into account, therefore double)
      // Default: 0.0 which means no adaptive sampling
      Double_t fToysInTails;
      // maximum no of toys
      // (taking weights into account, therefore double)
      Double_t fMaxToys;
      // tails
      Double_t fAdaptiveLowLimit;
      Double_t fAdaptiveHighLimit;

      const RooDataSet *fProtoData; // in dev
      
      ProofConfig *fProofConfig;   //!
      
      mutable NuisanceParametersSampler *fNuisanceParametersSampler; //!

      // objects below cache information and are mutable and non-persistent
      mutable RooArgSet* _allVars ; //! 
      mutable std::list<RooAbsPdf*> _pdfList ; //!
      mutable std::list<RooArgSet*> _obsList ; //!
      mutable std::list<RooAbsPdf::GenSpec*> _gsList ; //!      
      mutable RooAbsPdf::GenSpec* _gs1 ; //! GenSpec #1 
      mutable RooAbsPdf::GenSpec* _gs2 ; //! GenSpec #2
      mutable RooAbsPdf::GenSpec* _gs3 ; //! GenSpec #3
      mutable RooAbsPdf::GenSpec* _gs4 ; //! GenSpec #4
      
      static Bool_t fgAlwaysUseMultiGen ;  // Use PrepareMultiGen always
      Bool_t fUseMultiGen ; // Use PrepareMultiGen?

   protected:
   ClassDef(ToyMCSampler,3) // A simple implementation of the TestStatSampler interface
};
}


#endif
 ToyMCSampler.h:1
 ToyMCSampler.h:2
 ToyMCSampler.h:3
 ToyMCSampler.h:4
 ToyMCSampler.h:5
 ToyMCSampler.h:6
 ToyMCSampler.h:7
 ToyMCSampler.h:8
 ToyMCSampler.h:9
 ToyMCSampler.h:10
 ToyMCSampler.h:11
 ToyMCSampler.h:12
 ToyMCSampler.h:13
 ToyMCSampler.h:14
 ToyMCSampler.h:15
 ToyMCSampler.h:16
 ToyMCSampler.h:17
 ToyMCSampler.h:18
 ToyMCSampler.h:19
 ToyMCSampler.h:20
 ToyMCSampler.h:21
 ToyMCSampler.h:22
 ToyMCSampler.h:23
 ToyMCSampler.h:24
 ToyMCSampler.h:25
 ToyMCSampler.h:26
 ToyMCSampler.h:27
 ToyMCSampler.h:28
 ToyMCSampler.h:29
 ToyMCSampler.h:30
 ToyMCSampler.h:31
 ToyMCSampler.h:32
 ToyMCSampler.h:33
 ToyMCSampler.h:34
 ToyMCSampler.h:35
 ToyMCSampler.h:36
 ToyMCSampler.h:37
 ToyMCSampler.h:38
 ToyMCSampler.h:39
 ToyMCSampler.h:40
 ToyMCSampler.h:41
 ToyMCSampler.h:42
 ToyMCSampler.h:43
 ToyMCSampler.h:44
 ToyMCSampler.h:45
 ToyMCSampler.h:46
 ToyMCSampler.h:47
 ToyMCSampler.h:48
 ToyMCSampler.h:49
 ToyMCSampler.h:50
 ToyMCSampler.h:51
 ToyMCSampler.h:52
 ToyMCSampler.h:53
 ToyMCSampler.h:54
 ToyMCSampler.h:55
 ToyMCSampler.h:56
 ToyMCSampler.h:57
 ToyMCSampler.h:58
 ToyMCSampler.h:59
 ToyMCSampler.h:60
 ToyMCSampler.h:61
 ToyMCSampler.h:62
 ToyMCSampler.h:63
 ToyMCSampler.h:64
 ToyMCSampler.h:65
 ToyMCSampler.h:66
 ToyMCSampler.h:67
 ToyMCSampler.h:68
 ToyMCSampler.h:69
 ToyMCSampler.h:70
 ToyMCSampler.h:71
 ToyMCSampler.h:72
 ToyMCSampler.h:73
 ToyMCSampler.h:74
 ToyMCSampler.h:75
 ToyMCSampler.h:76
 ToyMCSampler.h:77
 ToyMCSampler.h:78
 ToyMCSampler.h:79
 ToyMCSampler.h:80
 ToyMCSampler.h:81
 ToyMCSampler.h:82
 ToyMCSampler.h:83
 ToyMCSampler.h:84
 ToyMCSampler.h:85
 ToyMCSampler.h:86
 ToyMCSampler.h:87
 ToyMCSampler.h:88
 ToyMCSampler.h:89
 ToyMCSampler.h:90
 ToyMCSampler.h:91
 ToyMCSampler.h:92
 ToyMCSampler.h:93
 ToyMCSampler.h:94
 ToyMCSampler.h:95
 ToyMCSampler.h:96
 ToyMCSampler.h:97
 ToyMCSampler.h:98
 ToyMCSampler.h:99
 ToyMCSampler.h:100
 ToyMCSampler.h:101
 ToyMCSampler.h:102
 ToyMCSampler.h:103
 ToyMCSampler.h:104
 ToyMCSampler.h:105
 ToyMCSampler.h:106
 ToyMCSampler.h:107
 ToyMCSampler.h:108
 ToyMCSampler.h:109
 ToyMCSampler.h:110
 ToyMCSampler.h:111
 ToyMCSampler.h:112
 ToyMCSampler.h:113
 ToyMCSampler.h:114
 ToyMCSampler.h:115
 ToyMCSampler.h:116
 ToyMCSampler.h:117
 ToyMCSampler.h:118
 ToyMCSampler.h:119
 ToyMCSampler.h:120
 ToyMCSampler.h:121
 ToyMCSampler.h:122
 ToyMCSampler.h:123
 ToyMCSampler.h:124
 ToyMCSampler.h:125
 ToyMCSampler.h:126
 ToyMCSampler.h:127
 ToyMCSampler.h:128
 ToyMCSampler.h:129
 ToyMCSampler.h:130
 ToyMCSampler.h:131
 ToyMCSampler.h:132
 ToyMCSampler.h:133
 ToyMCSampler.h:134
 ToyMCSampler.h:135
 ToyMCSampler.h:136
 ToyMCSampler.h:137
 ToyMCSampler.h:138
 ToyMCSampler.h:139
 ToyMCSampler.h:140
 ToyMCSampler.h:141
 ToyMCSampler.h:142
 ToyMCSampler.h:143
 ToyMCSampler.h:144
 ToyMCSampler.h:145
 ToyMCSampler.h:146
 ToyMCSampler.h:147
 ToyMCSampler.h:148
 ToyMCSampler.h:149
 ToyMCSampler.h:150
 ToyMCSampler.h:151
 ToyMCSampler.h:152
 ToyMCSampler.h:153
 ToyMCSampler.h:154
 ToyMCSampler.h:155
 ToyMCSampler.h:156
 ToyMCSampler.h:157
 ToyMCSampler.h:158
 ToyMCSampler.h:159
 ToyMCSampler.h:160
 ToyMCSampler.h:161
 ToyMCSampler.h:162
 ToyMCSampler.h:163
 ToyMCSampler.h:164
 ToyMCSampler.h:165
 ToyMCSampler.h:166
 ToyMCSampler.h:167
 ToyMCSampler.h:168
 ToyMCSampler.h:169
 ToyMCSampler.h:170
 ToyMCSampler.h:171
 ToyMCSampler.h:172
 ToyMCSampler.h:173
 ToyMCSampler.h:174
 ToyMCSampler.h:175
 ToyMCSampler.h:176
 ToyMCSampler.h:177
 ToyMCSampler.h:178
 ToyMCSampler.h:179
 ToyMCSampler.h:180
 ToyMCSampler.h:181
 ToyMCSampler.h:182
 ToyMCSampler.h:183
 ToyMCSampler.h:184
 ToyMCSampler.h:185
 ToyMCSampler.h:186
 ToyMCSampler.h:187
 ToyMCSampler.h:188
 ToyMCSampler.h:189
 ToyMCSampler.h:190
 ToyMCSampler.h:191
 ToyMCSampler.h:192
 ToyMCSampler.h:193
 ToyMCSampler.h:194
 ToyMCSampler.h:195
 ToyMCSampler.h:196
 ToyMCSampler.h:197
 ToyMCSampler.h:198
 ToyMCSampler.h:199
 ToyMCSampler.h:200
 ToyMCSampler.h:201
 ToyMCSampler.h:202
 ToyMCSampler.h:203
 ToyMCSampler.h:204
 ToyMCSampler.h:205
 ToyMCSampler.h:206
 ToyMCSampler.h:207
 ToyMCSampler.h:208
 ToyMCSampler.h:209
 ToyMCSampler.h:210
 ToyMCSampler.h:211
 ToyMCSampler.h:212
 ToyMCSampler.h:213
 ToyMCSampler.h:214
 ToyMCSampler.h:215
 ToyMCSampler.h:216
 ToyMCSampler.h:217
 ToyMCSampler.h:218
 ToyMCSampler.h:219
 ToyMCSampler.h:220
 ToyMCSampler.h:221
 ToyMCSampler.h:222
 ToyMCSampler.h:223
 ToyMCSampler.h:224
 ToyMCSampler.h:225
 ToyMCSampler.h:226
 ToyMCSampler.h:227
 ToyMCSampler.h:228
 ToyMCSampler.h:229
 ToyMCSampler.h:230
 ToyMCSampler.h:231
 ToyMCSampler.h:232
 ToyMCSampler.h:233
 ToyMCSampler.h:234
 ToyMCSampler.h:235
 ToyMCSampler.h:236
 ToyMCSampler.h:237
 ToyMCSampler.h:238
 ToyMCSampler.h:239
 ToyMCSampler.h:240
 ToyMCSampler.h:241
 ToyMCSampler.h:242
 ToyMCSampler.h:243
 ToyMCSampler.h:244
 ToyMCSampler.h:245
 ToyMCSampler.h:246
 ToyMCSampler.h:247
 ToyMCSampler.h:248
 ToyMCSampler.h:249
 ToyMCSampler.h:250
 ToyMCSampler.h:251
 ToyMCSampler.h:252
 ToyMCSampler.h:253
 ToyMCSampler.h:254
 ToyMCSampler.h:255
 ToyMCSampler.h:256
 ToyMCSampler.h:257
 ToyMCSampler.h:258
 ToyMCSampler.h:259
 ToyMCSampler.h:260
 ToyMCSampler.h:261
 ToyMCSampler.h:262
 ToyMCSampler.h:263
 ToyMCSampler.h:264
 ToyMCSampler.h:265
 ToyMCSampler.h:266
 ToyMCSampler.h:267
 ToyMCSampler.h:268
 ToyMCSampler.h:269
 ToyMCSampler.h:270
 ToyMCSampler.h:271
 ToyMCSampler.h:272
 ToyMCSampler.h:273
 ToyMCSampler.h:274
 ToyMCSampler.h:275
 ToyMCSampler.h:276
 ToyMCSampler.h:277
 ToyMCSampler.h:278
 ToyMCSampler.h:279
 ToyMCSampler.h:280
 ToyMCSampler.h:281
 ToyMCSampler.h:282
 ToyMCSampler.h:283
 ToyMCSampler.h:284
 ToyMCSampler.h:285
 ToyMCSampler.h:286
 ToyMCSampler.h:287
 ToyMCSampler.h:288
 ToyMCSampler.h:289
 ToyMCSampler.h:290
 ToyMCSampler.h:291
 ToyMCSampler.h:292
 ToyMCSampler.h:293
 ToyMCSampler.h:294
 ToyMCSampler.h:295
 ToyMCSampler.h:296
 ToyMCSampler.h:297
 ToyMCSampler.h:298
 ToyMCSampler.h:299
 ToyMCSampler.h:300
 ToyMCSampler.h:301
 ToyMCSampler.h:302
 ToyMCSampler.h:303
 ToyMCSampler.h:304
 ToyMCSampler.h:305
 ToyMCSampler.h:306
 ToyMCSampler.h:307
 ToyMCSampler.h:308
 ToyMCSampler.h:309
 ToyMCSampler.h:310
 ToyMCSampler.h:311
 ToyMCSampler.h:312
 ToyMCSampler.h:313
 ToyMCSampler.h:314
 ToyMCSampler.h:315
 ToyMCSampler.h:316
 ToyMCSampler.h:317
 ToyMCSampler.h:318
 ToyMCSampler.h:319
 ToyMCSampler.h:320
 ToyMCSampler.h:321
 ToyMCSampler.h:322
 ToyMCSampler.h:323
 ToyMCSampler.h:324
 ToyMCSampler.h:325
 ToyMCSampler.h:326
 ToyMCSampler.h:327