#include <Combine.h>

Public Member Functions
void	applyOptions (const boost::program_options::variables_map &vm)
	Combine ()
boost::program_options::options_description &	ioOptions ()
boost::program_options::options_description &	miscOptions ()
void	run (TString hlfFile, const std::string &dataset, double &limit, double &limitErr, int &iToy, TTree *tree, int nToys)
boost::program_options::options_description &	statOptions ()
Static Public Member Functions
static void	addBranch (const char name, void address, const char *leaflist)
	Add a branch to the output tree (for advanced use or debugging only)
static void	commitPoint (bool expected, float quantile)
	Save a point into the output tree. Usually if expected = false, quantile should be set to -1 (except e.g. for saveGrid option of HybridNew)
Private Member Functions
bool	mklimit (RooWorkspace w, RooStats::ModelConfig mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr)
Private Attributes
bool	compiledExpr_
float	expectSignal_
bool	generateBinnedWorkaround_
bool	guessGenMode_
bool	hintUsesStatOnly_
boost::program_options::options_description	ioOptions_
bool	makeTempDir_
float	mass_
boost::program_options::options_description	miscOptions_
std::string	modelConfigName_
std::string	modelConfigNameB_
bool	newGen_
bool	optSimPdf_
std::string	prior_
bool	rebuildSimPdf_
float	rMax_
float	rMin_
bool	saveToys_
bool	saveWorkspace_
boost::program_options::options_description	statOptions_
bool	toysFrequentist_
bool	toysNoSystematics_
bool	unbinned_
bool	validateModel_
std::string	workspaceName_
Static Private Attributes
static TTree *	tree_ = 0

Detailed Description

Definition at line 24 of file Combine.h.

Constructor & Destructor Documentation

Combine::Combine ( )

Definition at line 75 of file Combine.cc.

References cl, expectSignal_, ioOptions_, makeTempDir_, miscOptions_, modelConfigName_, modelConfigNameB_, newGen_, optSimPdf_, prior_, rebuildSimPdf_, rMax_, rMin_, statOptions_, withSystematics, and workspaceName_.

                 :
    statOptions_("Common statistics options"),
    ioOptions_("Common input-output options"),
    miscOptions_("Common miscellaneous options"),
    rMin_(std::numeric_limits<float>::quiet_NaN()), 
    rMax_(std::numeric_limits<float>::quiet_NaN()) {
    namespace po = boost::program_options;
    statOptions_.add_options() 
      ("systematics,S", po::value<bool>(&withSystematics)->default_value(true), "Add systematic uncertainties")
      ("cl,C",   po::value<float>(&cl)->default_value(0.95), "Confidence Level")
      ("rMin",   po::value<float>(&rMin_), "Override minimum value for signal strength")
      ("rMax",   po::value<float>(&rMax_), "Override maximum value for signal strength")
      ("prior",  po::value<std::string>(&prior_)->default_value("flat"), "Prior to use, for methods that require it and if it's not already in the input file: 'flat' (default), '1/sqrt(r)'")
      ("significance", "Compute significance instead of upper limit (works only for some methods)")
      ("lowerLimit",   "Compute the lower limit instead of the upper limit (works only for some methods)")
      ("hintStatOnly", "Ignore systematics when computing the hint")
      ("toysNoSystematics", "Generate all toys with the central value of the nuisance parameters, without fluctuating them")
      ("toysFrequentist",   "Generate all toys in the frequentist way. Note that this affects the toys generated by option '-t' that happen in the main loop, not the ones within the Hybrid(New) algorithm.")
      ("expectSignal", po::value<float>(&expectSignal_)->default_value(0.), "If set to non-zero, generate *signal* toys instead of background ones, with the specified signal strength.")
      ("unbinned,U", "Generate unbinned datasets instead of binned ones (works only for extended pdfs)")
      ("generateBinnedWorkaround", "Make binned datasets generating unbinned ones and then binnning them. Workaround for a bug in RooFit.")
      ;
    ioOptions_.add_options()
      ("saveWorkspace", "Save workspace to output root file")
      ("workspaceName,w", po::value<std::string>(&workspaceName_)->default_value("w"), "Workspace name, when reading it from or writing it to a rootfile.")
      ("modelConfigName",  po::value<std::string>(&modelConfigName_)->default_value("ModelConfig"), "ModelConfig name, when reading it from or writing it to a rootfile.")
      ("modelConfigNameB", po::value<std::string>(&modelConfigNameB_)->default_value("%s_bonly"), "Name of the ModelConfig for b-only hypothesis.\n"
                                                                                                  "If not present, it will be made from the singal model taking zero signal strength.\n"
                                                                                                  "A '%s' in the name will be replaced with the modelConfigName.")
      ("validateModel,V", "Perform some sanity checks on the model and abort if they fail.")
      ("saveToys",   "Save results of toy MC in output file")
      ;
    miscOptions_.add_options()
      ("newGenerator", po::value<bool>(&newGen_)->default_value(true), "Use new generator code for toys, fixes all issues with binned and mixed generation (equivalent of --newToyMC but affects the top-level toys from option '-t' instead of the ones within the HybridNew)")
      ("optimizeSimPdf", po::value<bool>(&optSimPdf_)->default_value(true), "Turn on special optimizations of RooSimultaneous. On by default, you can turn it off if it doesn't work for your workspace.")
      ("rebuildSimPdf", po::value<bool>(&rebuildSimPdf_)->default_value(false), "Rebuild simultaneous pdf from scratch to make sure constraints are correct (not needed in CMS workspaces)")
      ("compile", "Compile expressions instead of interpreting them")
      ("tempDir", po::value<bool>(&makeTempDir_)->default_value(false), "Run the program from a temporary directory (automatically on for text datacards or if 'compile' is activated)")
      ("guessGenMode", "Guess if to generate binned or unbinned based on dataset");
      ; 
}

Member Function Documentation

void Combine::addBranch	(	const char *	name,
		void *	address,
		const char *	leaflist
	)		`[static]`

Add a branch to the output tree (for advanced use or debugging only)

Definition at line 564 of file Combine.cc.

References tree_.

Referenced by MultiDimFit::initOnce(), and FitterAlgoBase::run().

                                                                             {
    tree_->Branch(name,address,leaflist);
}

void Combine::applyOptions ( const boost::program_options::variables_map & vm )

Definition at line 117 of file Combine.cc.

References compiledExpr_, gather_cfg::cout, doSignificance_, generateBinnedWorkaround_, guessGenMode_, hintUsesStatOnly_, lowerLimit_, makeTempDir_, mass_, modelConfigName_, modelConfigNameB_, saveToys_, saveWorkspace_, toysFrequentist_, toysNoSystematics_, unbinned_, validateModel_, and withSystematics.

                                                                      {
  if(withSystematics) {
    std::cout << ">>> including systematics" << std::endl;
  } else {
    std::cout << ">>> no systematics included" << std::endl;
  } 
  unbinned_ = vm.count("unbinned");
  generateBinnedWorkaround_ = vm.count("generateBinnedWorkaround");
  if (unbinned_ && generateBinnedWorkaround_) throw std::logic_error("You can't set generateBinnedWorkaround and unbinned options at the same time");
  guessGenMode_ = vm.count("guessGenMode");
  compiledExpr_ = vm.count("compile"); if (compiledExpr_) makeTempDir_ = true;
  doSignificance_ = vm.count("significance");
  lowerLimit_     = vm.count("lowerLimit");
  hintUsesStatOnly_ = vm.count("hintStatOnly");
  saveWorkspace_ = vm.count("saveWorkspace");
  toysNoSystematics_ = vm.count("toysNoSystematics");
  toysFrequentist_ = vm.count("toysFrequentist");
  if (toysNoSystematics_ && toysFrequentist_) throw std::logic_error("You can't set toysNoSystematics and toysFrequentist options at the same time");
  if (modelConfigNameB_.find("%s") != std::string::npos) {
      char modelBName[1024]; 
      sprintf(modelBName, modelConfigNameB_.c_str(), modelConfigName_.c_str());
      modelConfigNameB_ = modelBName;
  }
  mass_ = vm["mass"].as<float>();
  saveToys_ = vm.count("saveToys");
  validateModel_ = vm.count("validateModel");
}

void Combine::commitPoint	(	bool	expected,
		float	quantile
	)		`[static]`

Save a point into the output tree. Usually if expected = false, quantile should be set to -1 (except e.g. for saveGrid option of HybridNew)

Definition at line 557 of file Combine.cc.

References g_quantileExpected_, and tree_.

Referenced by MultiDimFit::doBox(), MultiDimFit::doContour2D(), MultiDimFit::doGrid(), MultiDimFit::doRandomPoints(), MultiDimFit::doSingles(), Asymptotic::getCLs(), AsymptoticNew::runLimit(), HybridNew::runLimit(), Asymptotic::runLimitExpected(), MultiDimFit::runSpecific(), and MaxLikelihoodFit::runSpecific().

                                                       {
    Float_t saveQuantile =  g_quantileExpected_;
    g_quantileExpected_ = quantile;
    tree_->Fill();
    g_quantileExpected_ = saveQuantile;
}

boost::program_options::options_description& Combine::ioOptions ( ) [inline]

Definition at line 29 of file Combine.h.

References ioOptions_.

{ return ioOptions_; }

boost::program_options::options_description& Combine::miscOptions ( ) [inline]

Definition at line 30 of file Combine.h.

References miscOptions_.

{ return miscOptions_; }

bool Combine::mklimit	(	RooWorkspace *	w,
		RooStats::ModelConfig *	mc_s,
		RooStats::ModelConfig *	mc_b,
		RooAbsData &	data,
		double &	limit,
		double &	limitErr
	)		`[private]`

Definition at line 145 of file Combine.cc.

References algo, dtNoiseDBValidation_cfg::cerr, exception, hintAlgo, hintUsesStatOnly_, run_regression::ret, LimitAlgo::run(), t_cpu_, t_real_, and withSystematics.

Referenced by run().

                                                                                                                                              {
  TStopwatch timer;
  bool ret = false;
  try {
    double hint = 0, hintErr = 0; bool hashint = false;
    if (hintAlgo) {
        if (hintUsesStatOnly_ && withSystematics) {
            withSystematics = false;
            hashint = hintAlgo->run(w, mc_s, mc_b, data, hint, hintErr, 0);
            withSystematics = true;
        } else {
            hashint = hintAlgo->run(w, mc_s, mc_b, data, hint, hintErr, 0);
        } 
    }
    limitErr = 0; // start with 0, as some algorithms don't compute it
    ret = algo->run(w, mc_s, mc_b, data, limit, limitErr, (hashint ? &hint : 0));    
  } catch (std::exception &ex) {
    std::cerr << "Caught exception " << ex.what() << std::endl;
    return false;
  }
  /* if ((ret == false) && (verbose > 3)) {
    std::cout << "Failed for method " << algo->name() << "\n";
    std::cout << "  --- DATA ---\n";
    utils::printRAD(&data);
    std::cout << "  --- MODEL ---\n";
    w->Print("V");
  } */
  timer.Stop(); t_cpu_ = timer.CpuTime()/60.; t_real_ = timer.RealTime()/60.;
  printf("Done in %.2f min (cpu), %.2f min (real)\n", t_cpu_, t_real_);
  return ret;
}

void Combine::run	(	TString	hlfFile,
		const std::string &	dataset,
		double &	limit,
		double &	limitErr,
		int &	iToy,
		TTree *	tree,
		int	nToys
	)

Definition at line 190 of file Combine.cc.

References algo, asimovutils::asimovDatasetWithFit(), dtNoiseDBValidation_cfg::cerr, utils::checkModel(), cl, compiledExpr_, gather_cfg::cout, dqm::qstatus::ERROR, expectSignal_, toymcoptutils::SimPdfGenInfo::generate(), toymcoptutils::SimPdfGenInfo::generateAsimov(), generateBinnedWorkaround_, utils::guessChannelMode(), guessGenMode_, instance, officialUncalib2006Production_0_4_0_cfg::isBinary, edm::detail::isnan(), MessageLogger_cff::limit, utils::makeNuisancePdf(), makeTempDir_, mass_, mklimit(), modelConfigName_, modelConfigNameB_, newGen_, AlCaHLTBitMon_ParallelJobs::options, optSimPdf_, download_sqlite_cfg::outputFile, funct::pow(), utils::printPdf(), utils::printRAD(), bookConverter::prior, prior_, alignCSCRings::pwd, readToysFromHere, utils::rebuildSimPdf(), rebuildSimPdf_, rMax_, rMin_, plotscripts::rms(), saveToys_, saveWorkspace_, utils::setAllConstant(), LimitAlgo::setNToys(), LimitAlgo::setToyNumber(), python::multivaluedict::sort(), mathSSE::sqrt(), ntuplemaker::status, toysFrequentist_, toysNoSystematics_, diffTreeTool::tree, tree_, unbinned_, validateModel_, w(), dqm::qstatus::WARNING, withSystematics, workspaceName_, and writeToysHere.

                                                                                                                               {
  ToCleanUp garbageCollect; // use this to close and delete temporary files

  TString tmpDir = "", tmpFile = "", pwd(gSystem->pwd());
  if (makeTempDir_) { 
      tmpDir = "roostats-XXXXXX"; tmpFile = "model";
      mkdtemp(const_cast<char *>(tmpDir.Data()));
      gSystem->cd(tmpDir.Data());
      garbageCollect.path = tmpDir.Data(); // request that we delete this dir when done
  } else if (!hlfFile.EndsWith(".hlf") && !hlfFile.EndsWith(".root")) {
      tmpFile = "roostats-XXXXXX";
      mktemp(const_cast<char *>(tmpFile.Data())); // somewhat unsafe, but I want to get a proper extension in the output file
  }

  bool isTextDatacard = false, isBinary = false;
  TString fileToLoad = (hlfFile[0] == '/' ? hlfFile : pwd+"/"+hlfFile);
  if (!boost::filesystem::exists(fileToLoad.Data())) throw std::invalid_argument(("File "+fileToLoad+" does not exist").Data());
  if (hlfFile.EndsWith(".hlf") ) {
    // nothing to do
  } else if (hlfFile.EndsWith(".root")) {
    isBinary = true;
  } else {
    TString txtFile = fileToLoad.Data();
    TString options = TString::Format(" -m %f -D %s", mass_, dataset.c_str());
    if (!withSystematics) options += " --stat ";
    if (compiledExpr_)    options += " --compiled ";
    if (verbose > 1)      options += TString::Format(" --verbose %d", verbose-1);
    //-- Text mode: old default
    //int status = gSystem->Exec("text2workspace.py "+options+" '"+txtFile+"' -o "+tmpFile+".hlf"); 
    //isTextDatacard = true; fileToLoad = tmpFile+".hlf";
    //-- Binary mode: new default 
    int status = gSystem->Exec("text2workspace.py "+options+" '"+txtFile+"' -b -o "+tmpFile+".root"); 
    isBinary = true; fileToLoad = tmpFile+".root";
    if (status != 0 || !boost::filesystem::exists(fileToLoad.Data())) {
        throw std::invalid_argument("Failed to convert the input datacard from LandS to RooStats format. The lines above probably contain more information about the error.");
    }
    garbageCollect.file = fileToLoad;
  }

  if (getenv("CMSSW_BASE")) {
      gSystem->AddIncludePath(TString::Format(" -I%s/src ", getenv("CMSSW_BASE")));
      if (verbose > 3) std::cout << "Adding " << getenv("CMSSW_BASE") << "/src to include path" << std::endl;
  }
  if (getenv("ROOFITSYS")) {
      gSystem->AddIncludePath(TString::Format(" -I%s/include ", getenv("ROOFITSYS")));
      if (verbose > 3) std::cout << "Adding " << getenv("ROOFITSYS") << "/include to include path" << std::endl;
  }

  if (verbose <= 2) RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  // Load the model, but going in a temporary directory to avoid polluting the current one with garbage from 'cexpr'
  RooWorkspace *w = 0; RooStats::ModelConfig *mc = 0, *mc_bonly = 0;
  std::auto_ptr<RooStats::HLFactory> hlf(0);
  if (isBinary) {
    TFile *fIn = TFile::Open(fileToLoad); 
    garbageCollect.tfile = fIn; // request that we close this file when done

    w = dynamic_cast<RooWorkspace *>(fIn->Get(workspaceName_.c_str()));
    if (w == 0) {  
        std::cerr << "Could not find workspace '" << workspaceName_ << "' in file " << fileToLoad << std::endl; fIn->ls(); 
        throw std::invalid_argument("Missing Workspace"); 
    }
    if (verbose > 3) { std::cout << "Input workspace '" << workspaceName_ << "': \n"; w->Print("V"); }
    RooRealVar *MH = w->var("MH");
    if (MH!=0) {
      if (verbose > 2) std::cerr << "Setting variable 'MH' in workspace to the higgs mass " << mass_ << std::endl;
      MH->setVal(mass_);
    }
    mc       = dynamic_cast<RooStats::ModelConfig *>(w->genobj(modelConfigName_.c_str()));
    mc_bonly = dynamic_cast<RooStats::ModelConfig *>(w->genobj(modelConfigNameB_.c_str()));
    if (mc == 0) {  
        std::cerr << "Could not find ModelConfig '" << modelConfigName_ << "' in workspace '" << workspaceName_ << "' in file " << fileToLoad << std::endl;
        throw std::invalid_argument("Missing ModelConfig"); 
    } else if (verbose > 3) { std::cout << "Workspace has a ModelConfig for signal called '" << modelConfigName_ << "', with contents:\n"; mc->Print("V"); }
    if (verbose > 3) { std::cout << "Input ModelConfig '" << modelConfigName_ << "': \n"; mc->Print("V"); }
    const RooArgSet *POI = mc->GetParametersOfInterest();
    if (POI == 0 || POI->getSize() == 0) throw std::invalid_argument("ModelConfig '"+modelConfigName_+"' does not define parameters of interest.");
    if (POI->getSize() > 1) std::cerr << "ModelConfig '" << modelConfigName_ << "' defines more than one parameter of interest. This is not supported in some statistical methods." << std::endl;
    if (mc->GetObservables() == 0) throw std::invalid_argument("ModelConfig '"+modelConfigName_+"' does not define observables.");
    if (mc->GetPdf() == 0) throw std::invalid_argument("ModelConfig '"+modelConfigName_+"' does not define a pdf.");
    if (rebuildSimPdf_ && typeid(*mc->GetPdf()) == typeid(RooSimultaneous)) {
        RooSimultaneous *newpdf = utils::rebuildSimPdf(*mc->GetObservables(), dynamic_cast<RooSimultaneous*>(mc->GetPdf()));
        w->import(*newpdf);
        mc->SetPdf(*newpdf);
    }
    if (optSimPdf_ && typeid(*mc->GetPdf()) == typeid(RooSimultaneous)) {
        RooSimultaneousOpt *optpdf = new RooSimultaneousOpt(static_cast<RooSimultaneous&>(*mc->GetPdf()), TString(mc->GetPdf()->GetName())+"_opt");
        w->import(*optpdf);
        mc->SetPdf(*optpdf);
    }
    if (mc_bonly == 0) {
        std::cerr << "Missing background ModelConfig '" << modelConfigNameB_ << "' in workspace '" << workspaceName_ << "' in file " << fileToLoad << std::endl;
        std::cerr << "Will make one from the signal ModelConfig '" << modelConfigName_ << "' setting signal strenth '" << POI->first()->GetName() << "' to zero"  << std::endl;
        w->factory("_zero_[0]");
        RooCustomizer make_model_s(*mc->GetPdf(),"_model_bonly_");
        make_model_s.replaceArg(*POI->first(), *w->var("_zero_"));
        RooAbsPdf *model_b = dynamic_cast<RooAbsPdf *>(make_model_s.build()); 
        model_b->SetName("_model_bonly_");
        w->import(*model_b);
        mc_bonly = new RooStats::ModelConfig(*mc);
        mc_bonly->SetPdf(*model_b);
    }
  } else {
    hlf.reset(new RooStats::HLFactory("factory", fileToLoad));
    w = hlf->GetWs();
    if (w == 0) {
        std::cerr << "Could not read HLF from file " <<  (hlfFile[0] == '/' ? hlfFile : pwd+"/"+hlfFile) << std::endl;
        return;
    }
    RooRealVar *MH = w->var("MH");
    if (MH==0) {
      std::cerr << "Could not find MH var in workspace '" << workspaceName_ << "' in file " << fileToLoad << std::endl;
      throw std::invalid_argument("Missing MH"); 
    }
    MH->setVal(mass_);
    if (w->set("observables") == 0) throw std::invalid_argument("The model must define a RooArgSet 'observables'");
    if (w->set("POI")         == 0) throw std::invalid_argument("The model must define a RooArgSet 'POI' for the parameters of interest");
    if (w->pdf("model_b")     == 0) throw std::invalid_argument("The model must define a RooAbsPdf 'model_b'");
    if (w->pdf("model_s")     == 0) throw std::invalid_argument("The model must define a RooAbsPdf 'model_s'");

    // create ModelConfig
    mc = new RooStats::ModelConfig(modelConfigName_.c_str(),"signal",w);
    mc->SetPdf(*w->pdf("model_s"));
    mc->SetObservables(*w->set("observables"));
    mc->SetParametersOfInterest(*w->set("POI"));
    if (w->set("nuisances"))         mc->SetNuisanceParameters(*w->set("nuisances"));
    if (w->set("globalObservables")) mc->SetGlobalObservables(*w->set("globalObservables"));
    if (w->pdf("prior")) mc->SetNuisanceParameters(*w->pdf("prior"));
    w->import(*mc, modelConfigName_.c_str());

    mc_bonly = new RooStats::ModelConfig(modelConfigNameB_.c_str(),"background",w);
    mc_bonly->SetPdf(*w->pdf("model_b"));
    mc_bonly->SetObservables(*w->set("observables"));
    mc_bonly->SetParametersOfInterest(*w->set("POI"));
    if (w->set("nuisances"))         mc_bonly->SetNuisanceParameters(*w->set("nuisances"));
    if (w->set("globalObservables")) mc_bonly->SetGlobalObservables(*w->set("globalObservables"));
    if (w->pdf("prior")) mc_bonly->SetNuisanceParameters(*w->pdf("prior"));
    w->import(*mc_bonly, modelConfigNameB_.c_str());
  }
  gSystem->cd(pwd);

  if (verbose <= 2) RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);

  const RooArgSet * observables = mc->GetObservables();     // not null
  const RooArgSet * POI = mc->GetParametersOfInterest();     // not null
  const RooArgSet * nuisances = mc->GetNuisanceParameters(); // note: may be null
  if (dynamic_cast<RooRealVar*>(POI->first()) == 0) throw std::invalid_argument("First parameter of interest is not a RooRealVar");

  if (w->data(dataset.c_str()) == 0) {
    if (isTextDatacard) { // that's ok: the observables are pre-set to the observed values
      RooDataSet *data_obs = new RooDataSet(dataset.c_str(), dataset.c_str(), *observables); 
      data_obs->add(*observables);
      w->import(*data_obs);
    } else {
      std::cout << "Dataset " << dataset.c_str() << " not found." << std::endl;
    }
  }

  if (verbose < -1) {
      RooMsgService::instance().setStreamStatus(0,kFALSE);
      RooMsgService::instance().setStreamStatus(1,kFALSE);
      RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
  }


  if (!isnan(rMin_)) ((RooRealVar*)POI->first())->setMin(rMin_);
  if (!isnan(rMax_)) ((RooRealVar*)POI->first())->setMax(rMax_);

  if (mc->GetPriorPdf() == 0) {
      if (prior_ == "flat") {
          RooAbsPdf *prior = new RooUniform("prior","prior",*POI);
          w->import(*prior);
          mc->SetPriorPdf(*prior);
      } else if (prior_ == "1/sqrt(r)") {
          w->factory(TString::Format("EXPR::prior(\\\"1/sqrt(@0)\\\",%s)", POI->first()->GetName()));
          mc->SetPriorPdf(*w->pdf("prior"));
      } else if (!prior_.empty() && w->pdf(prior_.c_str()) != 0) {
          std::cout << "Will use prior '" << prior_ << "' in from the input workspace" << std::endl;
          mc->SetPriorPdf(*w->pdf(prior_.c_str()));
      } else {
          std::cerr << "Unknown prior '" << prior_ << "'. It's not 'flat' '1/sqrt(r)' or the name of a pdf in the model.\n" << std::endl;
          throw std::invalid_argument("Bad prior");
      }
  }

  if (withSystematics && nuisances == 0) {
      std::cout << "The signal model has no nuisance parameters. Please run the limit tool with no systematics (option -S 0)." << std::endl;
      std::cout << "To make things easier, I will assume you have done it." << std::endl;
      withSystematics = false;
  } else if (!withSystematics && nuisances != 0) {
    std::cout << "Will set nuisance parameters to constants: " ;
    utils::setAllConstant(*nuisances, true);
  }

  bool validModel = validateModel_ ? utils::checkModel(*mc, false) : true;
  if (validateModel_ && verbose) std::cout << "Sanity checks on the model: " << (validModel ? "OK" : "FAIL") << std::endl;

  // make sure these things are set consistently with what we expect
  if (mc->GetNuisanceParameters() && withSystematics) utils::setAllConstant(*mc->GetNuisanceParameters(), false);
  if (mc->GetGlobalObservables()) utils::setAllConstant(*mc->GetGlobalObservables(), true);
  if (mc_bonly->GetNuisanceParameters() && withSystematics) utils::setAllConstant(*mc_bonly->GetNuisanceParameters(), false);
  if (mc_bonly->GetGlobalObservables()) utils::setAllConstant(*mc_bonly->GetGlobalObservables(), true);


  w->saveSnapshot("clean", w->allVars());

  if (saveWorkspace_) {
    w->SetName(workspaceName_.c_str());
    outputFile->WriteTObject(w,workspaceName_.c_str());
  }

  tree_ = tree;

  bool isExtended = mc_bonly->GetPdf()->canBeExtended();
  RooAbsData *dobs = w->data(dataset.c_str());
  RooAbsPdf  *genPdf = expectSignal_ > 0 ? mc->GetPdf() : mc_bonly->GetPdf();
  toymcoptutils::SimPdfGenInfo newToyMC(*genPdf, *observables, !unbinned_); RooRealVar *weightVar_ = 0;

  if (guessGenMode_ && genPdf->InheritsFrom("RooSimultaneous") && (dobs != 0)) {
      utils::guessChannelMode(dynamic_cast<RooSimultaneous&>(*mc->GetPdf()), *dobs, verbose);
      utils::guessChannelMode(dynamic_cast<RooSimultaneous&>(*mc_bonly->GetPdf()), *dobs, 0);
  }
  if (expectSignal_ > 0) { 
      ((RooRealVar*)POI->first())->setVal(expectSignal_); 
  }
  if (nToys <= 0) { // observed or asimov
    iToy = nToys;
    if (iToy == -1) {   
        if (newGen_) {
            if (toysFrequentist_) {
                w->saveSnapshot("reallyClean", w->allVars());
                if (dobs == 0) throw std::invalid_argument("Frequentist Asimov datasets can't be generated without a real dataset to fit");
                RooArgSet gobsAsimov;
                dobs = asimovutils::asimovDatasetWithFit(mc, *dobs, gobsAsimov, expectSignal_, verbose);
                if (mc->GetGlobalObservables()) {
                    RooArgSet gobs(*mc->GetGlobalObservables());
                    gobs = gobsAsimov;
                    utils::setAllConstant(*mc->GetParametersOfInterest(), false);
                    w->saveSnapshot("clean", w->allVars());
                }
            } else {
                dobs = newToyMC.generateAsimov(weightVar_); // as simple as that
            }
        } else if (isExtended) {
            if (unbinned_) {
                throw std::invalid_argument("Asimov datasets can only be generated binned");
            } else {
                dobs = genPdf->generateBinned(*observables,RooFit::Extended(),RooFit::Asimov());
            }
        } else {
          dobs = genPdf->generate(*observables,1,RooFit::Asimov());
        }
    } else if (dobs == 0) {
      std::cerr << "No observed data '" << dataset << "' in the workspace. Cannot compute limit.\n" << std::endl;
      return;
    }
    if (saveToys_) {
        writeToysHere->WriteTObject(dobs, TString::Format("toy_asimov%g", expectSignal_));
      }
    std::cout << "Computing limit starting from " << (iToy == 0 ? "observation" : "expected outcome") << std::endl;
    if (verbose > (isExtended ? 3 : 2)) utils::printRAD(dobs);
    if (mklimit(w,mc,mc_bonly,*dobs,limit,limitErr)) tree->Fill();
  }
  
  
  std::vector<double> limitHistory;
  std::auto_ptr<RooAbsPdf> nuisancePdf;
  if (nToys > 0) {
    double expLimit = 0;
    unsigned int nLimits = 0;
    w->loadSnapshot("clean");
    RooDataSet *systDs = 0;
    if (withSystematics && !toysNoSystematics_ && (readToysFromHere == 0)) {
      if (nuisances == 0) throw std::logic_error("Running with systematics enabled, but nuisances not defined.");
      nuisancePdf.reset(utils::makeNuisancePdf(expectSignal_ ? *mc : *mc_bonly));
      if (toysFrequentist_) {
          if (mc->GetGlobalObservables() == 0) throw std::logic_error("Cannot use toysFrequentist with no global observables");
          w->saveSnapshot("reallyClean", w->allVars());
          {
              if (dobs == 0) throw std::logic_error("Cannot use toysFrequentist with no input dataset");
              CloseCoutSentry sentry(verbose < 3);
              genPdf->fitTo(*dobs, RooFit::Save(1), RooFit::Minimizer("Minuit2","minimize"), RooFit::Strategy(0), RooFit::Hesse(0), RooFit::Constrain(*(expectSignal_ ?mc:mc_bonly)->GetNuisanceParameters()));
          }
          w->saveSnapshot("clean", w->allVars());
          systDs = nuisancePdf->generate(*mc->GetGlobalObservables(), nToys);
      } else {
          systDs = nuisancePdf->generate(*nuisances, nToys);
      } 
    }
    std::auto_ptr<RooArgSet> vars(genPdf->getVariables());
    algo->setNToys(nToys);
    for (iToy = 1; iToy <= nToys; ++iToy) {
      algo->setToyNumber(iToy-1);
      RooAbsData *absdata_toy = 0;
      if (readToysFromHere == 0) {
        w->loadSnapshot("clean");
        if (verbose > 3) utils::printPdf(*mc_bonly);
        if (withSystematics && !toysNoSystematics_) {
          *vars = *systDs->get(iToy-1);
          if (toysFrequentist_) w->saveSnapshot("clean", w->allVars());
          if (verbose > 3) utils::printPdf(*mc_bonly);
        }
        if (expectSignal_) ((RooRealVar*)POI->first())->setVal(expectSignal_);
        std::cout << "Generate toy " << iToy << "/" << nToys << std::endl;
        if (isExtended) {
          if (newGen_) {
              absdata_toy = newToyMC.generate(weightVar_); // as simple as that
          } else if (unbinned_) {
              absdata_toy = genPdf->generate(*observables,RooFit::Extended());
          } else if (generateBinnedWorkaround_) {
              std::auto_ptr<RooDataSet> unbinn(genPdf->generate(*observables,RooFit::Extended()));
              absdata_toy = new RooDataHist("toy","binned toy", *observables, *unbinn);
          } else {
              absdata_toy = genPdf->generateBinned(*observables,RooFit::Extended());
          }
        } else {
          RooDataSet *data_toy = genPdf->generate(*observables,1);
          absdata_toy = data_toy;
        }
      } else {
        absdata_toy = dynamic_cast<RooAbsData *>(readToysFromHere->Get(TString::Format("toys/toy_%d",iToy)));
        if (absdata_toy == 0) {
          std::cerr << "Toy toy_"<<iToy<<" not found in " << readToysFromHere->GetName() << ". List follows:\n";
          readToysFromHere->ls();
          return;
        }
      }
      if (verbose > (isExtended ? 3 : 2)) utils::printRAD(absdata_toy);
      w->loadSnapshot("clean");
      //if (verbose > 1) utils::printPdf(w, "model_b");
      if (mklimit(w,mc,mc_bonly,*absdata_toy,limit,limitErr)) {
        tree->Fill();
        ++nLimits;
        expLimit += limit; 
        limitHistory.push_back(limit);
      }
      if (saveToys_) {
        writeToysHere->WriteTObject(absdata_toy, TString::Format("toy_%d", iToy));
      }
      delete absdata_toy;
    }
    if (weightVar_) delete weightVar_;
    expLimit /= nLimits;
    double rms = 0;
    for (std::vector<double>::const_iterator itl = limitHistory.begin(); itl != limitHistory.end(); ++itl) {
        rms += pow(*itl-expLimit, 2);
    }
    if (nLimits > 1) {
        rms = sqrt(rms/(nLimits-1)/nLimits);
        cout << "mean   expected limit: r < " << expLimit << " +/- " << rms << " @ " << cl*100 << "%CL (" <<nLimits << " toyMC)" << endl;
    } else {
        cout << "mean   expected limit: r < " << expLimit << " @ " << cl*100 << "%CL (" <<nLimits << " toyMC)" << endl;
    }
    sort(limitHistory.begin(), limitHistory.end());
    if (nLimits > 0) {
        double medianLimit = (nLimits % 2 == 0 ? 0.5*(limitHistory[nLimits/2-1]+limitHistory[nLimits/2]) : limitHistory[nLimits/2]);
        cout << "median expected limit: r < " << medianLimit << " @ " << cl*100 << "%CL (" <<nLimits << " toyMC)" << endl;
        double hi68 = limitHistory[min<int>(nLimits-1,  ceil(0.84  * nLimits))];
        double lo68 = limitHistory[min<int>(nLimits-1, floor(0.16  * nLimits))];
        double hi95 = limitHistory[min<int>(nLimits-1,  ceil(0.975 * nLimits))];
        double lo95 = limitHistory[min<int>(nLimits-1, floor(0.025 * nLimits))];
        cout << "   68% expected band : " << lo68 << " < r < " << hi68 << endl;
        cout << "   95% expected band : " << lo95 << " < r < " << hi95 << endl;
    }
  }

}