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>
 
ClassImp(RooSimSplitGenContext);
 
 
////////////////////////////////////////////////////////////////////////////////
/// 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() ;
  _fracThresh = new double[_numPdf+1] ;
  _fracThresh[0] = 0 ;
 
  // Generate index category and all registered PDFS
  _allVarsPdf.add(allPdfVars) ;
  Int_t i(1) ;
  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);
 
    // Fill fraction threshold array
    _fracThresh[i] = _fracThresh[i-1] + pdf->expectedEvents(&allPdfVars) ;
    i++ ;
  }
 
  for(i=0 ; i<_numPdf ; i++) {
    _fracThresh[i] /= _fracThresh[_numPdf] ;
  }
 
  // 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()
{
  delete[] _fracThresh ;
  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) ;
  if (extendedMode ) {
    Int_t i(0) ;
    for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
      RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
      //nGen[i] = Int_t(pdf->expectedEvents(&_allVarsPdf)+0.5) ;
      nGen[i] = pdf->expectedEvents(&_allVarsPdf) ;
      i++ ;
    }
 
  } else {
    Int_t i(1) ;
    _fracThresh[0] = 0 ;
    for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
      RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
      _fracThresh[i] = _fracThresh[i-1] + pdf->expectedEvents(&_allVarsPdf) ;
      i++ ;
    }
    for(i=0 ; i<_numPdf ; i++) {
      _fracThresh[i] /= _fracThresh[_numPdf] ;
    }
 
    // Determine from that total number of events to be generated for each component
    double nGenSoFar(0) ;
    while (nGenSoFar<nEvents) {
      double rand = RooRandom::uniform() ;
      i=0 ;
      for (i=0 ; i<_numPdf ; i++) {
   if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
     nGen[i]++ ;
     nGenSoFar++ ;
     break ;
   }
      }
    }
  }
 
 
 
  // 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)
{
  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);
}