RooSimSplitGenContext.cxx

Go to the documentation of this file.

/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * 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)             *
 *****************************************************************************/
 
/**
\file RooSimSplitGenContext.cxx
\class RooSimSplitGenContext
\ingroup Roofitcore
 
Efficient implementation of the generator context
specific for RooSimultaneous PDFs when generating more than one of the
component pdfs.
**/
 
#include "RooSimSplitGenContext.h"
#include "RooSimultaneous.h"
#include "RooRealProxy.h"
#include "RooDataSet.h"
#include "Roo1DTable.h"
#include "RooCategory.h"
#include "RooMsgService.h"
#include "RooRandom.h"
#include "RooGlobalFunc.h"
 
#include <iostream>
#include <string>
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor of specialized generator context for RooSimultaneous p.d.f.s. This
/// context creates a dedicated context for each component p.d.f.s and delegates
/// generation of events to the appropriate component generator context
 
RooSimSplitGenContext::RooSimSplitGenContext(const RooSimultaneous &model, const RooArgSet &vars, bool verbose, bool autoBinned, const char* binnedTag) :
  RooAbsGenContext(model,vars,nullptr,nullptr,verbose), _pdf(&model)
{
  // Determine if we are requested to generate the index category
  RooAbsCategoryLValue const& idxCat = model.indexCat();
  RooArgSet pdfVars(vars) ;
 
  RooArgSet allPdfVars(pdfVars) ;
 
  RooArgSet catsAmongAllVars;
  allPdfVars.selectCommon(model.flattenedCatList(), catsAmongAllVars);
 
  if(catsAmongAllVars.size() != model.flattenedCatList().size()) {
      oocoutE(_pdf,Generation) << "RooSimSplitGenContext::ctor(" << GetName() << ") ERROR: This context must"
                << " generate all components of the index category" << std::endl ;
      _isValid = false ;
      _numPdf = 0 ;
      // coverity[UNINIT_CTOR]
      return ;
  }
 
  // We must extended likelihood to determine the relative fractions of the components
  _idxCatName = idxCat.GetName() ;
  if (!model.canBeExtended()) {
    oocoutE(_pdf,Generation) << "RooSimSplitGenContext::RooSimSplitGenContext(" << GetName() << "): All components of the simultaneous PDF "
              << "must be extended PDFs. Otherwise, it is impossible to calculate the number of events to be generated per component." << std::endl ;
    _isValid = false ;
    _numPdf = 0 ;
    // coverity[UNINIT_CTOR]
    return ;
  }
 
  // Initialize fraction threshold array (used only in extended mode)
  _numPdf = model._pdfProxyList.GetSize() ;
 
  // Generate index category and all registered PDFS
  _allVarsPdf.add(allPdfVars) ;
  for(auto * proxy : static_range_cast<RooRealProxy*>(model._pdfProxyList)) {
    auto pdf = static_cast<RooAbsPdf*>(proxy->absArg());
 
    // Create generator context for this PDF
    std::unique_ptr<RooArgSet> compVars{pdf->getObservables(pdfVars)};
    RooAbsGenContext* cx = pdf->autoGenContext(*compVars,nullptr,nullptr,verbose,autoBinned,binnedTag) ;
 
    const auto state = idxCat.lookupIndex(proxy->name());
 
    cx->SetName(proxy->name()) ;
    _gcList.push_back(cx) ;
    _gcIndex.push_back(state);
  }
 
  // Clone the index category
  if(RooArgSet(model.indexCat()).snapshot(_idxCatSet, true)) {
    oocoutE(_pdf,Generation) << "RooSimSplitGenContext::RooSimSplitGenContext(" << GetName() << ") Couldn't deep-clone index category, abort," << std::endl ;
    throw std::string("RooSimSplitGenContext::RooSimSplitGenContext() Couldn't deep-clone index category, abort") ;
  }
  _idxCat = static_cast<RooAbsCategoryLValue*>(_idxCatSet.find(model.indexCat().GetName()));
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Destructor. Delete all owned subgenerator contexts
 
RooSimSplitGenContext::~RooSimSplitGenContext()
{
  for (RooAbsGenContext *item : _gcList) {
    delete item;
  }
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Attach the index category clone to the given event buffer
 
void RooSimSplitGenContext::attach(const RooArgSet& args)
{
  if (_idxCat->isDerived()) {
    _idxCat->recursiveRedirectServers(args) ;
  }
 
  // Forward initGenerator call to all components
  for (RooAbsGenContext *item : _gcList) {
    item->attach(args) ;
  }
 
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Perform one-time initialization of generator context
 
void RooSimSplitGenContext::initGenerator(const RooArgSet &theEvent)
{
  // Attach the index category clone to the event
  if (_idxCat->isDerived()) {
    _idxCat->recursiveRedirectServers(theEvent) ;
  } else {
    _idxCat = static_cast<RooAbsCategoryLValue*>(theEvent.find(_idxCat->GetName())) ;
  }
 
  // Forward initGenerator call to all components
  for (RooAbsGenContext *item : _gcList) {
    item->initGenerator(theEvent) ;
  }
 
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
 
RooDataSet* RooSimSplitGenContext::generate(double nEvents, bool skipInit, bool extendedMode)
{
  if(!isValid()) {
    coutE(Generation) << ClassName() << "::" << GetName() << ": context is not valid" << std::endl;
    return nullptr;
  }
 
 
  // Calculate the expected number of events if necessary
  if(nEvents <= 0) {
    nEvents= _expectedEvents;
  }
  coutI(Generation) << ClassName() << "::" << GetName() << ":generate: will generate "
          << nEvents << " events" << std::endl;
 
  if (_verbose) Print("v") ;
 
  // Perform any subclass implementation-specific initialization
  // Can be skipped if this is a rerun with an identical configuration
  if (!skipInit) {
    initGenerator(_theEvent);
  }
 
  // Generate lookup table from expected event counts
  std::vector<double> nGen(_numPdf);
  std::vector<double> nExpected;
  nExpected.reserve(_numPdf) ;
  double nExpectedTotal = 0.;
  for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
    RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
    nExpected.push_back(pdf->expectedEvents(&_allVarsPdf));
    nExpectedTotal += nExpected.back();
  }
 
  // We don't randomize events in two cases:
  //  1. When the generation is extended, each component pdf will already
  //     randomize the expected number of events so we don't need to do it here.
  //  2. If we want to create an expected Asimov dataset.
  if (extendedMode || _expectedData) {
    nGen = nExpected;
  } else {
    // Determine from that total number of events to be generated for each component
    double nGenSoFar(0) ;
    while (nGenSoFar<nEvents) {
      double rand = RooRandom::uniform() * nExpectedTotal;
      double cumsum = 0;
      for (int i=0 ; i<_numPdf ; i++) {
        if (rand >= cumsum && rand < cumsum + nExpected[i]) {
          nGen[i]++ ;
          nGenSoFar++ ;
          break ;
        }
        cumsum += nExpected[i];
      }
    }
  }
 
 
 
  // Now loop over states
  std::map<std::string,RooAbsData*> dataMap ;
  Int_t icomp(0) ;
  for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
 
    // Calculate number of events to generate for this state
    if (_gcList[icomp]) {
      dataMap[proxy->GetName()] = _gcList[icomp]->generate(nGen[icomp],skipInit,extendedMode) ;
    }
 
    icomp++ ;
  }
 
  // Put all datasets together in a composite-store RooDataSet that links and owns the component datasets
  RooDataSet* hmaster = new RooDataSet("hmaster","hmaster",_allVarsPdf,RooFit::Index(static_cast<RooCategory&>(*_idxCat)),RooFit::Link(dataMap),RooFit::OwnLinked()) ;
  return hmaster ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Forward to components
 
void RooSimSplitGenContext::setExpectedData(bool flag)
{
  _expectedData = flag;
  for(RooAbsGenContext *elem : _gcList) {
    elem->setExpectedData(flag) ;
  }
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// this method is empty because it is not used by this context
 
RooDataSet* RooSimSplitGenContext::createDataSet(const char* , const char* , const RooArgSet& )
{
  return nullptr;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// this method is empty because it is not used in this type of context
 
void RooSimSplitGenContext::generateEvent(RooArgSet &, Int_t )
{
  assert(0) ;
}
 
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// this method is empty because proto datasets are not supported by this context
 
void RooSimSplitGenContext::setProtoDataOrder(Int_t* )
{
  assert(0) ;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Detailed printing interface
 
void RooSimSplitGenContext::printMultiline(std::ostream &os, Int_t content, bool verbose, TString indent) const
{
  RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
  os << indent << "--- RooSimSplitGenContext ---" << std::endl ;
  os << indent << "Using PDF ";
  _pdf->printStream(os,kName|kArgs|kClassName,kSingleLine,indent);
}