    /*****************************************************************************
     * Project: RooFit                                                           *
     * Package: RooFitCore                                                       *
     *    File: $Id: RooMCStudy.cc,v 1.30 2005/06/20 15:44:54 wverkerke Exp $
     * Authors:                                                                  *
     *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
     *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
     *                                                                           *
     * Copyright (c) 2000-2005, Regents of the University of California          *
     *                          and Stanford University. All rights reserved.    *
     *                                                                           *
     * Redistribution and use in source and binary forms,                        *
     * with or without modification, are permitted according to the terms        *
     * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
     *****************************************************************************/

    // -- CLASS DESCRIPTION [MISC] --
    // RooMCStudy is a help class to facilitate Monte Carlo studies
    // such as 'goodness-of-fit' studies, that involve fitting a PDF
    // to multiple toy Monte Carlo sets generated from the same PDF
    // or another PDF.
    //
    // Given a fit PDF and a generator PDF, RooMCStudy can produce
    // large numbers of toyMC samples and/or fit these samples
    // and acculumate the final parameters of each fit in a dataset.
    //
    // Additional plotting routines simplify the task of plotting
    // the distribution of the minimized likelihood, each parameters fitted value,
    // fitted error and pull distribution.

    #include "RooFit.h"

    #include "RooMCStudy.h"
    #include "RooMCStudy.h"

    #include "RooGenContext.h"
    #include "RooAbsPdf.h"
    #include "RooDataSet.h"
    #include "RooDataHist.h"
    #include "RooRealVar.h"
    #include "RooFitResult.h"
    #include "RooErrorVar.h"
    #include "RooFormulaVar.h"
    #include "RooArgList.h"
    #include "RooPlot.h"
    #include "RooGenericPdf.h"
    #include "RooRandom.h"
    #include "RooCmdConfig.h"
    #include "RooGlobalFunc.h"
    #include "RooPullVar.h"

    ClassImp(RooMCStudy)
      ;

    RooMCStudy::RooMCStudy(const RooAbsPdf& model, const RooArgSet& observables,
       		       RooCmdArg arg1, RooCmdArg arg2,
       		       RooCmdArg arg3,RooCmdArg arg4,RooCmdArg arg5,
       		       RooCmdArg arg6,RooCmdArg arg7,RooCmdArg arg8)

      // Construct Monte Carlo Study Manager. This class automates generating data from a given PDF,
      // fitting the PDF to that data and accumulating the fit statistics.
      //
      // The constructor accepts the following arguments
      //
      // model       -- The PDF to be studied
      // observables -- The variables of the PDF to be considered the observables
      //
      // FitModel(const RooAbsPdf&)        -- The PDF for fitting, if it is different from the PDF for generating
      // ConditionalObservables
      //           (const RooArgSet& set)  -- The set of observables that the PDF should _not_ be normalized over
      // Binned(Bool_t flag)               -- Bin the dataset before fitting it. Speeds up fitting of large data samples
      // FitOptions(....)                  -- Options to be used for fitting. All named arguments inside FitOptions()
      //                                                   are passed to RooAbsPdf::fitTo();
      // Verbose(Bool_t flag)              -- Activate informational messages in event generation phase
      // Extended(Bool_t flag)             -- Determine number of events for each sample anew from a Poisson distribution
      // ProtoData(const RooDataSet&,
      //                 Bool_t randOrder) -- Prototype data for the event generation. If the randOrder flag is
      //                                      set, the order of the dataset will be re-randomized for each generation
      //                                      cycle to protect against systematic biases if the number of generated
      //                                      events does not exactly match the number of events in the prototype dataset
      //                                      at the cost of reduced precision
      //                                      with mu equal to the specified number of events
    {
      // Stuff all arguments in a list
      RooLinkedList cmdList;
      cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;

      // Select the pdf-specific commands
      RooCmdConfig pc(Form("RooMCStudy::RooMCStudy(%s)",model.GetName())) ;

      pc.defineObject("fitModel","FitModel",0,0) ;
      pc.defineObject("condObs","ProjectedDependents",0,0) ;
      pc.defineObject("protoData","PrototypeData",0,0) ;
      pc.defineInt("randProtoData","PrototypeData",0,0) ;
      pc.defineInt("verboseGen","Verbose",0,0) ;
      pc.defineInt("extendedGen","Extended",0,0) ;
      pc.defineInt("binGenData","Binned",0,0) ;
      pc.defineString("fitOpts","FitOptions",0,"") ;
      pc.defineInt("dummy","FitOptArgs",0,0) ;
      pc.defineMutex("FitOptions","FitOptArgs") ;

      // Process and check varargs
      pc.process(cmdList) ;
      if (!pc.ok(kTRUE)) {
        // WVE do something here
        return ;
      }

      // Save fit command options
      if (pc.hasProcessed("FitOptArgs")) {
        RooCmdArg* fitOptArg = static_cast<RooCmdArg*>(cmdList.FindObject("FitOptArgs")) ;
        for (Int_t i=0 ; i<fitOptArg->subArgs().GetSize() ;i++) {
          _fitOptList.Add(new RooCmdArg(static_cast<RooCmdArg&>(*fitOptArg->subArgs().At(i)))) ;
        }
      }

      // Decode command line arguments
      _verboseGen = pc.getInt("verboseGen") ;
      _extendedGen = pc.getInt("extendedGen") ;
      _binGenData = pc.getInt("binGenData") ;
      _randProto = pc.getInt("randProtoData") ;

      _genModel = const_cast<RooAbsPdf*>(&model) ;
      RooAbsPdf* fitModel = static_cast<RooAbsPdf*>(pc.getObject("fitModel",0)) ;
      _fitModel = fitModel ? fitModel : _genModel ;

      _genProtoData = static_cast<RooDataSet*>(pc.getObject("protoData",0)) ;
      if (pc.getObject("condObs",0)) {
        _projDeps.add(static_cast<RooArgSet&>(*pc.getObject("condObs",0))) ;
      }

      _dependents.add(observables) ;

      _allDependents.add(_dependents) ;
      _fitOptions = pc.getString("fitOpts") ;
      _canAddFitResults = kTRUE ;

      if (_extendedGen && _genProtoData && !_randProto) {
        cout << "RooMCStudy::RooMCStudy: WARNING Using generator option 'e' (Poisson distribution of #events) together " << endl
    	 << "                        with a prototype dataset implies incomplete sampling or oversampling of proto data." << endl
    	 << "                        Use option \"r\" to randomize prototype dataset order and thus to randomize" << endl
    	 << "                        the set of over/undersampled prototype events for each generation cycle." << endl ;
      }

      _genContext = _genModel->genContext(_dependents,_genProtoData,0,_verboseGen) ;
      RooArgSet* tmp = _genModel->getParameters(&_dependents) ;
      _genParams = (RooArgSet*) tmp->snapshot(kFALSE) ;
      delete tmp ;

      // Store list of parameters and save initial values separately
      _fitParams = _fitModel->getParameters(&_dependents) ;
      _fitInitParams = (RooArgSet*) _fitParams->snapshot(kTRUE) ;

      _nExpGen = _extendedGen ? _genModel->expectedEvents(&_dependents) : 0 ;

      // Place holder for NLL
      _nllVar = new RooRealVar("NLL","-log(Likelihood)",0) ;

      // Create data set containing parameter values, errors and pulls
      RooArgSet tmp2(*_fitParams) ;
      tmp2.add(*_nllVar) ;

      // Mark all variable to store their errors in the dataset
      tmp2.setAttribAll("StoreError",kTRUE) ;
      tmp2.setAttribAll("StoreAsymError",kTRUE) ;
      _fitParData = new RooDataSet("fitParData","Fit Parameters DataSet",tmp2) ;
      tmp2.setAttribAll("StoreError",kFALSE) ;
      tmp2.setAttribAll("StoreAsymError",kFALSE) ;

      // Append proto variables to allDependents
      if (_genProtoData) {
        _allDependents.add(*_genProtoData->get(),kTRUE) ;
      }

    }

     RooMCStudy::RooMCStudy(const RooAbsPdf& genModel, const RooAbsPdf& fitModel,
       		       const RooArgSet& dependents, const char* genOptions,
       		       const char* fitOptions, const RooDataSet* genProtoData,
       		       const RooArgSet& projDeps) :
      _genModel((RooAbsPdf*)&genModel),
      _genProtoData(genProtoData),
      _projDeps(projDeps),
      _dependents(dependents),
      _allDependents(dependents),
      _fitModel((RooAbsPdf*)&fitModel),
      _fitOptions(fitOptions),
      _canAddFitResults(kTRUE)
    {
      // Constructor with a generator and fit model. Both models may point
      // to the same object. The 'dependents' set of variables is generated
      // in the generator phase. The optional prototype dataset is passed to
      // the generator
      //
      // Available generator options
      //  v  - Verbose
      //  e  - Extended: use Poisson distribution for Nevts generated
      //
      // Available fit options
      //  See RooAbsPdf::fitTo()
      //

      // Decode generator options
      TString genOpt(genOptions) ;
      genOpt.ToLower() ;
      _verboseGen = genOpt.Contains("v") ;
      _extendedGen = genOpt.Contains("e") ;
      _binGenData = genOpt.Contains("b") ;
      _randProto = genOpt.Contains("r") ;

      if (_extendedGen && genProtoData && !_randProto) {
        cout << "RooMCStudy::RooMCStudy: WARNING Using generator option 'e' (Poisson distribution of #events) together " << endl
       	 << "                        with a prototype dataset implies incomplete sampling or oversampling of proto data." << endl
       	 << "                        Use option \"r\" to randomize prototype dataset order and thus to randomize" << endl
    	 << "                        the set of over/undersampled prototype events for each generation cycle." << endl ;
      }

      _genContext = genModel.genContext(dependents,genProtoData,0,_verboseGen) ;
      RooArgSet* tmp = genModel.getParameters(&dependents) ;
      _genParams = (RooArgSet*) tmp->snapshot(kFALSE) ;
      delete tmp ;

      // Store list of parameters and save initial values separately
      _fitParams = fitModel.getParameters(&dependents) ;
      _fitInitParams = (RooArgSet*) _fitParams->snapshot(kTRUE) ;

      _nExpGen = _extendedGen ? genModel.expectedEvents(&dependents) : 0 ;

      // Place holder for NLL
      _nllVar = new RooRealVar("NLL","-log(Likelihood)",0) ;

      // Create data set containing parameter values, errors and pulls
      RooArgSet tmp2(*_fitParams) ;
      tmp2.add(*_nllVar) ;

      // Mark all variable to store their errors in the dataset
      tmp2.setAttribAll("StoreError",kTRUE) ;
      tmp2.setAttribAll("StoreAsymError",kTRUE) ;
      _fitParData = new RooDataSet("fitParData","Fit Parameters DataSet",tmp2) ;
      tmp2.setAttribAll("StoreError",kFALSE) ;
      tmp2.setAttribAll("StoreAsymError",kFALSE) ;

      // Append proto variables to allDependents
      if (genProtoData) {
        _allDependents.add(*genProtoData->get(),kTRUE) ;
      }
    }

     RooMCStudy::~RooMCStudy()
    {
      // Destructor

      _genDataList.Delete() ;
      _fitResList.Delete() ;
      _fitOptList.Delete() ;
      delete _fitParData ;
      delete _fitParams ;
      delete _genParams ;
      delete _genContext ;
      delete _nllVar ;
    }

     Bool_t RooMCStudy::run(Bool_t generate, Bool_t fit, Int_t nSamples, Int_t nEvtPerSample, Bool_t keepGenData, const char* asciiFilePat)
    {
      // Run engine. Generate and/or fit, according to flags, 'nSamples' samples of 'nEvtPerSample' events.
      // If keepGenData is set, all generated data sets will be kept in memory and can be accessed
      // later via genData().
      //
      // When generating, data sets will be written out in ascii form if the pattern string is supplied
      // The pattern, which is a template for sprintf, should look something like "data/toymc_%04d.dat"
      // and should contain one integer field that encodes the sample serial number.
      //
      // When fitting only, data sets may optionally be read from ascii files, using the same file
      // pattern.
      //

      while(nSamples--) {

        cout << "RooMCStudy::run: " ;
        if (generate) cout << "Generating " ;
        if (generate && fit) cout << "and " ;
        if (fit) cout << "fitting " ;
        cout << "sample " << nSamples << endl ;

        RooDataSet* genSample = 0;
        if (generate) {
          // Generate sample
          Int_t nEvt(nEvtPerSample) ;

          if (_extendedGen) {
    	nEvt = RooRandom::randomGenerator()->Poisson(nEvtPerSample==0?_nExpGen:nEvtPerSample) ;
          }

          if (_randProto && _genProtoData && _genProtoData->numEntries()!=nEvt) {
    	cout << "RooMCStudy: (Re)randomizing event order in prototype dataset (Nevt=" << nEvt << ")" << endl ;
              Int_t* newOrder = _genModel->randomizeProtoOrder(_genProtoData->numEntries(),nEvt) ;
              _genContext->setProtoDataOrder(newOrder) ;
              delete[] newOrder ;
          }

          genSample = _genContext->generate(nEvt) ;

        } else if (asciiFilePat && &asciiFilePat) {
          // Load sample from ASCII file
          char asciiFile[1024] ;
          sprintf(asciiFile,asciiFilePat,nSamples) ;
          RooArgList depList(_allDependents) ;
          genSample = RooDataSet::read(asciiFile,depList,"q") ;

        } else {

          // Load sample from internal list
          genSample = (RooDataSet*) _genDataList.At(nSamples) ;
          if (!genSample) {
       	cout << "RooMCStudy::run: WARNING: Sample #" << nSamples << " not loaded, skipping" << endl ;
       	continue ;
          }
        }

        if (fit) fitSample(genSample) ;

        // Optionally write to ascii file
        if (generate && asciiFilePat && *asciiFilePat) {
          char asciiFile[1024] ;
          sprintf(asciiFile,asciiFilePat,nSamples) ;
          genSample->write(asciiFile) ;
        }

        // Add to list or delete
        if (keepGenData) {
          _genDataList.Add(genSample) ;
        } else {
          delete genSample ;
        }
      }

      _canAddFitResults = kFALSE ;
      if (fit) calcPulls() ;
      return kFALSE ;
    }

     Bool_t RooMCStudy::generateAndFit(Int_t nSamples, Int_t nEvtPerSample, Bool_t keepGenData, const char* asciiFilePat)
    {
      // Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
      // If keepGenData is set, all generated data sets will be kept in memory and can be accessed
      // later via genData().
      //
      // Data sets will be written out is ascii form if the pattern string is supplied.
      // The pattern, which is a template for sprintf, should look something like "data/toymc_%04d.dat"
      // and should contain one integer field that encodes the sample serial number.
      //

      // Clear any previous data in memory
      _fitResList.Delete() ;
      _genDataList.Delete() ;
      _fitParData->reset() ;

      return run(kTRUE,kTRUE,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
    }

     Bool_t RooMCStudy::generate(Int_t nSamples, Int_t nEvtPerSample, Bool_t keepGenData, const char* asciiFilePat)
    {
      // Generate 'nSamples' samples of 'nEvtPerSample' events.
      // If keepGenData is set, all generated data sets will be kept in memory
      // and can be accessed later via genData().
      //
      // Data sets will be written out in ascii form if the pattern string is supplied.
      // The pattern, which is a template for sprintf, should look something like "data/toymc_%04d.dat"
      // and should contain one integer field that encodes the sample serial number.
      //

      // Clear any previous data in memory
      _genDataList.Delete() ;

      return run(kTRUE,kFALSE,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
    }

     Bool_t RooMCStudy::fit(Int_t nSamples, const char* asciiFilePat)
    {
      // Fit 'nSamples' datasets, which are read from ASCII files.
      //
      // The ascii file pattern, which is a template for sprintf, should look something like "data/toymc_%04d.dat"
      // and should contain one integer field that encodes the sample serial number.
      //

      // Clear any previous data in memory
      _fitResList.Delete() ;
      _fitParData->reset() ;

      return run(kFALSE,kTRUE,nSamples,0,kFALSE,asciiFilePat) ;
    }

     Bool_t RooMCStudy::fit(Int_t nSamples, TList& dataSetList)
    {
      // Fit 'nSamples' datasets, as supplied in 'dataSetList'
      //

      // Clear any previous data in memory
      _fitResList.Delete() ;
      _genDataList.Delete() ;
      _fitParData->reset() ;

      // Load list of data sets
      TIterator* iter = dataSetList.MakeIterator() ;
      RooAbsData* gset ;
      while((gset=(RooAbsData*)iter->Next())) {
        _genDataList.Add(gset) ;
      }
      delete iter ;

      return run(kFALSE,kTRUE,nSamples,0,kTRUE,0) ;
    }

     Bool_t RooMCStudy::fitSample(RooAbsData* genSample)
    {
      // Fit given dataset with fit model. If fit
      // converges (TMinuit status code zero)
      // The fit results are appended to the fit results
      // dataset
      //
      // If the fit option "r" is supplied, the RooFitResult
      // objects will always be saved, regardless of the
      // fit status. RooFitResults objects can be retrieved
      // later via fitResult().
      //

      // Reset all fit parameters to their initial values
      *_fitParams = *_fitInitParams ;

      // Fit model to data set
      TString fitOpt2(_fitOptions) ; fitOpt2.Append("r") ;

      // Optionally bin dataset before fitting
      RooAbsData* data ;
      if (_binGenData) {
        RooArgSet* depList = _fitModel->getObservables(genSample) ;
        data = new RooDataHist(genSample->GetName(),genSample->GetTitle(),*depList,*genSample) ;
        delete depList ;
      } else {
        data = genSample ;
      }

      RooFitResult* fr ;
      if (_fitOptList.GetSize()==0) {
        if (_projDeps.getSize()>0) {
          fr = (RooFitResult*) _fitModel->fitTo(*data,_projDeps,fitOpt2) ;
        } else {
          fr = (RooFitResult*) _fitModel->fitTo(*data,fitOpt2) ;
        }
      } else {
        RooCmdArg save  = RooFit::Save() ;
        RooCmdArg condo = RooFit::ConditionalObservables(_projDeps) ;
        RooLinkedList fitOptList(_fitOptList) ;
        fitOptList.Add(&save) ;
        if (_projDeps.getSize()>0) {
          fitOptList.Add(&condo) ;
        }
        fr = (RooFitResult*) _fitModel->fitTo(*data,fitOptList) ;
      }

      // If fit converged, store parameters and NLL
      Bool_t ok = (fr->status()==0) ;
      if (ok) {
        _nllVar->setVal(fr->minNll()) ;
        RooArgSet tmp(*_fitParams) ;
        tmp.add(*_nllVar) ;
        _fitParData->add(tmp) ;
      }

      // Store fit result if requested by user
      Bool_t userSaveRequest = kFALSE ;
      if (_fitOptList.GetSize()>0) {
        if (_fitOptList.FindObject("Save")) userSaveRequest = kTRUE ;
      } else {
        if (_fitOptions.Contains("r")) userSaveRequest = kTRUE ;
      }

      if (userSaveRequest) {
        _fitResList.Add(fr) ;
      } else {
        delete fr ;
      }

      if (_binGenData) delete data ;

      return !ok ;
    }

     Bool_t RooMCStudy::addFitResult(const RooFitResult& fr)
    {
      if (!_canAddFitResults) {
        cout << "RooMCStudy::addFitResult: ERROR cannot add fit results in current state" << endl ;
        return kTRUE ;
      }

      // Transfer contents of fit result to fitParams ;
      *_fitParams = RooArgSet(fr.floatParsFinal()) ;

      // If fit converged, store parameters and NLL
      Bool_t ok = (fr.status()==0) ;
      if (ok) {
        _nllVar->setVal(fr.minNll()) ;
        RooArgSet tmp(*_fitParams) ;
        tmp.add(*_nllVar) ;
        _fitParData->add(tmp) ;
      }

      // Store fit result if requested by user
      if (_fitOptions.Contains("r")) {
        _fitResList.Add((TObject*)&fr) ;
      }

      return kFALSE ;
    }

     void RooMCStudy::calcPulls()
    {
      // Calculate the pulls for all fit parameters in
      // the fit results data set, and add them to that dataset

      TIterator* iter = _fitParams->createIterator()  ;
      RooRealVar* par ;
      while((par=(RooRealVar*)iter->Next())) {

        RooErrorVar* err = par->errorVar() ;
        _fitParData->addColumn(*err) ;

        TString name(par->GetName()), title(par->GetTitle()) ;
        name.Append("pull") ;
        title.Append(" Pull") ;
        RooAbsReal* genParOrig = (RooAbsReal*)_genParams->find(par->GetName()) ;
        if (genParOrig) {
          RooAbsReal* genPar = (RooAbsReal*) genParOrig->Clone("truth") ;
          RooPullVar pull(name,title,*par,*genPar) ;

          _fitParData->addColumn(pull) ;
          delete genPar ;
        }
      }
      delete iter ;

    }

     const RooDataSet& RooMCStudy::fitParDataSet()
    {
      // Return the fit parameter dataset
      if (_canAddFitResults) {
        calcPulls() ;
        _canAddFitResults = kFALSE ;
      }

      return *_fitParData ;
    }

     const RooArgSet* RooMCStudy::fitParams(Int_t sampleNum) const
    {
      // Return an argset with the fit parameters for the given sample number
      // NB: The fit parameters are only stored for successfull fits,
      //     thus the maximum sampleNum can be less that the number
      //     of generated samples and if so, the indeces will
      //     be out of synch with genData() and fitResult()

      // Check if sampleNum is in range
      if (sampleNum<0 || sampleNum>=_fitParData->numEntries()) {
        cout << "RooMCStudy::fitParams: ERROR, invalid sample number: " << sampleNum << endl ;
        return 0 ;
      }

      return _fitParData->get(sampleNum) ;
    }

     const RooFitResult* RooMCStudy::fitResult(Int_t sampleNum) const
    {
      // Return the fit result object of the fit to given sample

      // Check if sampleNum is in range
      if (sampleNum<0 || sampleNum>=_fitResList.GetSize()) {
        cout << "RooMCStudy::fitResult: ERROR, invalid sample number: " << sampleNum << endl ;
        return 0 ;
      }

      // Retrieve fit result object
      const RooFitResult* fr = (RooFitResult*) _fitResList.At(sampleNum) ;
      if (fr) {
        return fr ;
      } else {
        cout << "RooMCStudy::fitResult: ERROR, no fit result saved for sample "
       	 << sampleNum << ", did you use the 'r; fit option?" << endl ;
      }
      return 0 ;
    }

     const RooDataSet* RooMCStudy::genData(Int_t sampleNum) const
    {
      // Return the given generated dataset

      // Check that generated data was saved
      if (_genDataList.GetSize()==0) {
        cout << "RooMCStudy::genData() ERROR, generated data was not saved" << endl ;
        return 0 ;
      }

      // Check if sampleNum is in range
      if (sampleNum<0 || sampleNum>=_genDataList.GetSize()) {
        cout << "RooMCStudy::genData() ERROR, invalid sample number: " << sampleNum << endl ;
        return 0 ;
      }

      return (RooDataSet*) _genDataList.At(sampleNum) ;
    }

     RooPlot* RooMCStudy::plotParamOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
       				 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
    {
      // Plot the distribution of the fitted value of the given parameter on the specified frame
      // Any specified named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options

      _fitParData->plotOn(frame,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
      return frame ;
    }

     RooPlot* RooMCStudy::plotParam(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
       			       const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
    {
      // Plot the distribution of the fitted value of the given parameter on a newly created frame.
      //
      // This function accepts the following optional arguments
      // FrameRange(double lo, double hi) -- Set range of frame to given specification
      // FrameBins(int bins)              -- Set default number of bins of frame to given number
      // Frame(...)                       -- Pass supplied named arguments to RooAbsRealLValue::frame() function. See frame() function
      //                                     for list of allowed arguments
      //
      // If no frame specifications are given, the AutoRange() feature will be used to set the range
      // Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options

      // Stuff all arguments in a list
      RooLinkedList cmdList;
      cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;

      RooPlot* frame = makeFrameAndPlotCmd(param, cmdList) ;
      if (frame) {
        _fitParData->plotOn(frame, cmdList) ;
      }

      return frame ;
    }

     RooPlot* RooMCStudy::plotNLL(const RooCmdArg& arg1, const RooCmdArg& arg2,
                         const RooCmdArg& arg3, const RooCmdArg& arg4,
                         const RooCmdArg& arg5, const RooCmdArg& arg6,
                         const RooCmdArg& arg7, const RooCmdArg& arg8)
    {
      // Plot the distribution of the -log(l) values on a newly created frame.
      //
      // This function accepts the following optional arguments
      // FrameRange(double lo, double hi) -- Set range of frame to given specification
      // FrameBins(int bins)              -- Set default number of bins of frame to given number
      // Frame(...)                       -- Pass supplied named arguments to RooAbsRealLValue::frame() function. See frame() function
      //                                     for list of allowed arguments
      //
      // If no frame specifications are given, the AutoRange() feature will be used to set the range
      // Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options

      return plotParam(*_nllVar,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
    }

     RooPlot* RooMCStudy::plotError(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
                         const RooCmdArg& arg3, const RooCmdArg& arg4,
                         const RooCmdArg& arg5, const RooCmdArg& arg6,
                         const RooCmdArg& arg7, const RooCmdArg& arg8)
    {
      // Plot the distribution of the fit errors for the specified parameter on a newly created frame.
      //
      // This function accepts the following optional arguments
      // FrameRange(double lo, double hi) -- Set range of frame to given specification
      // FrameBins(int bins)              -- Set default number of bins of frame to given number
      // Frame(...)                       -- Pass supplied named arguments to RooAbsRealLValue::frame() function. See frame() function
      //                                     for list of allowed arguments
      //
      // If no frame specifications are given, the AutoRange() feature will be used to set the range
      // Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options

      if (_canAddFitResults) {
        calcPulls() ;
        _canAddFitResults=kFALSE ;
      }

      RooErrorVar* evar = param.errorVar() ;
      RooRealVar* evar_rrv = static_cast<RooRealVar*>(evar->createFundamental()) ;
      RooPlot* frame = plotParam(*evar_rrv,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
      delete evar_rrv ;
      delete evar ;
      return frame ;
    }

     RooPlot* RooMCStudy::plotPull(const RooRealVar& param, const RooCmdArg& arg1, const RooCmdArg& arg2,
                         const RooCmdArg& arg3, const RooCmdArg& arg4,
                         const RooCmdArg& arg5, const RooCmdArg& arg6,
                         const RooCmdArg& arg7, const RooCmdArg& arg8)
    {
      // Plot the distribution of pull values for the specified parameter on a newly created frame. If asymmetric
      // errors are calculated in the fit (by MINOS) those will be used in the pull calculation
      //
      // This function accepts the following optional arguments
      // FrameRange(double lo, double hi) -- Set range of frame to given specification
      // FrameBins(int bins)              -- Set default number of bins of frame to given number
      // Frame(...)                       -- Pass supplied named arguments to RooAbsRealLValue::frame() function. See frame() function
      //                                     for list of allowed arguments
      // FitGauss(Bool_t flag)            -- Add a gaussian fit to the frame
      //
      // If no frame specifications are given, the AutoSymRange() feature will be used to set the range
      // Any other named argument is passed to the RooAbsData::plotOn() call. See that function for allowed options

      // Stuff all arguments in a list
      RooLinkedList cmdList;
      cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
      cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ;  cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;

      TString name(param.GetName()), title(param.GetTitle()) ;
      name.Append("pull") ; title.Append(" Pull") ;
      RooRealVar pvar(name,title,-100,100) ;
      pvar.setBins(100) ;

      RooPlot* frame = makeFrameAndPlotCmd(pvar, cmdList, kTRUE) ;
      if (frame) {

        // Pick up optonal FitGauss command from list
        RooCmdConfig pc(Form("RooMCStudy::plotPull(%s)",_genModel->GetName())) ;
        pc.defineInt("fitGauss","FitGauss",0,0) ;
        pc.allowUndefined() ;
        pc.process(cmdList) ;
        Bool_t fitGauss=pc.getInt("fitGauss") ;

        // Pass stripped command list to plotOn()
        pc.stripCmdList(cmdList,"FitGauss") ;
        _fitParData->plotOn(frame,cmdList) ;

        // Add Gaussian fit if requested
        if (fitGauss) {
          RooRealVar pullMean("pullMean","Mean of pull",0,-100,100) ;
          RooRealVar pullSigma("pullSigma","Width of pull",1,0,5) ;
          RooGenericPdf pullGauss("pullGauss","Gaussian of pull",
    			      "exp(-0.5*(@0-@1)*(@0-@1)/(@2*@2))",
    			      RooArgSet(pvar,pullMean,pullSigma)) ;
          pullGauss.fitTo(*_fitParData,RooFit::Minos(0),RooFit::PrintLevel(-1)) ;
          pullGauss.plotOn(frame) ;
          pullGauss.paramOn(frame,_fitParData) ;
        }
      }
      return frame ; ;
    }

     RooPlot* RooMCStudy::makeFrameAndPlotCmd(const RooRealVar& param, RooLinkedList& cmdList, Bool_t symRange) const
    {

      // Select the frame-specific commands
      RooCmdConfig pc(Form("RooMCStudy::plotParam(%s)",_genModel->GetName())) ;
      pc.defineInt("nbins","FrameBins",0,0) ;
      pc.defineDouble("xlo","FrameRange",0,0) ;
      pc.defineDouble("xhi","FrameRange",1,0) ;
      pc.defineInt("dummy","FrameArgs",0,0) ;
      pc.defineMutex("FrameBins","FrameArgs") ;
      pc.defineMutex("FrameRange","FrameArgs") ;

      // Process and check varargs
      pc.allowUndefined() ;
      pc.process(cmdList) ;
      if (!pc.ok(kTRUE)) {
        return 0 ;
      }

      // Make frame according to specs
      Int_t nbins = pc.getInt("nbins") ;
      Double_t xlo = pc.getDouble("xlo") ;
      Double_t xhi = pc.getDouble("xhi") ;
      RooPlot* frame ;

      if (pc.hasProcessed("FrameArgs")) {
        // Explicit frame arguments are given, pass them on
        RooCmdArg* frameArg = static_cast<RooCmdArg*>(cmdList.FindObject("FrameArgs")) ;
        frame = param.frame(frameArg->subArgs()) ;
      } else {
        // FrameBins, FrameRange or none are given, build custom frame command list
        RooCmdArg bins = RooFit::Bins(nbins) ;
        RooCmdArg range = RooFit::Range(xlo,xhi) ;
        RooCmdArg autor = symRange ? RooFit::AutoSymRange(*_fitParData,0.2) : RooFit::AutoRange(*_fitParData,0.2) ;
        RooLinkedList frameCmdList ;

        if (pc.hasProcessed("FrameBins")) frameCmdList.Add(&bins) ;
        if (pc.hasProcessed("FrameRange")) {
          frameCmdList.Add(&range) ;
        } else {
          frameCmdList.Add(&autor) ;
        }
        frame = param.frame(frameCmdList) ;
      }

      // Filter frame command from list and pass on to plotOn()
      pc.stripCmdList(cmdList,"FrameBins,FrameRange,FrameArgs") ;

      return frame ;
    }

     RooPlot* RooMCStudy::plotNLL(Double_t lo, Double_t hi, Int_t nBins)
    {
      // Create a RooPlot of the NLL distribution in the range lo-hi
      // with 'nBins' bins

      RooPlot* frame = _nllVar->frame(lo,hi,nBins) ;

      _fitParData->plotOn(frame) ;
      return frame ;
    }

     RooPlot* RooMCStudy::plotError(const RooRealVar& param, Double_t lo, Double_t hi, Int_t nbins)
    {
      // Create a RooPlot of the distribution of the fitted errors of the given parameter.
      // The range lo-hi is plotted in nbins bins
      if (_canAddFitResults) {
        calcPulls() ;
        _canAddFitResults=kFALSE ;
      }

      RooErrorVar* evar = param.errorVar() ;
      RooPlot* frame = evar->frame(lo,hi,nbins) ;
      _fitParData->plotOn(frame) ;

      delete evar ;
      return frame ;
    }

     RooPlot* RooMCStudy::plotPull(const RooRealVar& param, Double_t lo, Double_t hi, Int_t nbins, Bool_t fitGauss)
    {
      // Create a RooPlot of the pull distribution for the given parameter.
      // The range lo-hi is plotted in nbins.
      // If fitGauss is set, an unbinned max. likelihood fit of the distribution to a Gaussian model
      // is performed. The fit result is overlaid on the returned RooPlot and a box with the fitted
      // mean and sigma is added.

      if (_canAddFitResults) {
        calcPulls() ;
        _canAddFitResults=kFALSE ;
      }

      TString name(param.GetName()), title(param.GetTitle()) ;
      name.Append("pull") ; title.Append(" Pull") ;
      RooRealVar pvar(name,title,lo,hi) ;
      pvar.setBins(nbins) ;

      RooPlot* frame = pvar.frame() ;
      _fitParData->plotOn(frame) ;

      if (fitGauss) {
        RooRealVar pullMean("pullMean","Mean of pull",0,lo,hi) ;
        RooRealVar pullSigma("pullSigma","Width of pull",1,0,5) ;
        RooGenericPdf pullGauss("pullGauss","Gaussian of pull",
    			    "exp(-0.5*(@0-@1)*(@0-@1)/(@2*@2))",
    			    RooArgSet(pvar,pullMean,pullSigma)) ;
        pullGauss.fitTo(*_fitParData,"mh") ;
        pullGauss.plotOn(frame) ;
        pullGauss.paramOn(frame,_fitParData) ;
      }

      return frame ;
    }