RooAddGenContext.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 RooAddGenContext.cxx
\class RooAddGenContext
\ingroup Roofitcore
 
Efficient implementation of the
generator context specific for RooAddPdf PDFs. The strategy
of RooAddGenContext is to defer generation of each component
to a dedicated generator context for that component and to
randomly choose one of those context to generate an event,
with a probability proportional to its associated coefficient.
**/
 
#include "RooAddGenContext.h"
 
#include "RooDataSet.h"
#include "RooRandom.h"
 
#include <sstream>
 
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor
 
RooAddGenContext::RooAddGenContext(const RooAddPdf &model, const RooArgSet &vars,
               const RooDataSet *prototype, const RooArgSet* auxProto,
               bool verbose) :
  RooAbsGenContext(model,vars,prototype,auxProto,verbose), _isModel(false)
{
  cxcoutI(Generation) << "RooAddGenContext::ctor() setting up event special generator context for sum p.d.f. " << model.GetName()
         << " for generation of observable(s) " << vars ;
  if (prototype) ccxcoutI(Generation) << " with prototype data for " << *prototype->get() ;
  if (auxProto && !auxProto->empty())  ccxcoutI(Generation) << " with auxiliary prototypes " << *auxProto ;
  ccxcoutI(Generation) << std::endl;
 
  // Constructor. Build an array of generator contexts for each product component PDF
  _pdfSet = std::make_unique<RooArgSet>();
  RooArgSet(model).snapshot(*_pdfSet, true);
  _pdf = static_cast<RooAddPdf*>(_pdfSet->find(model.GetName())) ;
 
  // Fix normalization set of this RooAddPdf
  if (prototype)
    {
      RooArgSet coefNSet(vars) ;
      coefNSet.add(*prototype->get()) ;
      _pdf->fixAddCoefNormalization(coefNSet,false) ;
    }
 
  _nComp = model._pdfList.size() ;
  _coefThresh.resize(_nComp+1);
  _vars = std::make_unique<RooArgSet>();
  vars.snapshot(*_vars, false);
 
  for (const auto arg : model._pdfList) {
    auto pdf = dynamic_cast<const RooAbsPdf *>(arg);
    if (!pdf) {
      coutF(Generation) << "Cannot generate events from an object that is not a PDF.\n\t"
          << "The offending object is a " << arg->ClassName() << " named '" << arg->GetName() << "'." << std::endl;
      throw std::invalid_argument("Trying to generate events from on object that is not a PDF.");
    }
 
    _gcList.emplace_back(pdf->genContext(vars,prototype,auxProto,verbose));
  }
 
  (static_cast<RooAddPdf*>(_pdf))->getProjCache(_vars.get()) ;
  _pdf->recursiveRedirectServers(_theEvent) ;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor
 
RooAddGenContext::RooAddGenContext(const RooAddModel &model, const RooArgSet &vars,
               const RooDataSet *prototype, const RooArgSet* auxProto,
               bool verbose) :
  RooAbsGenContext(model,vars,prototype,auxProto,verbose), _isModel(true)
{
  cxcoutI(Generation) << "RooAddGenContext::ctor() setting up event special generator context for sum resolution model " << model.GetName()
         << " for generation of observable(s) " << vars ;
  if (prototype) ccxcoutI(Generation) << " with prototype data for " << *prototype->get() ;
  if (auxProto && !auxProto->empty())  ccxcoutI(Generation) << " with auxiliary prototypes " << *auxProto ;
  ccxcoutI(Generation) << std::endl;
 
  // Constructor. Build an array of generator contexts for each product component PDF
  _pdfSet = std::make_unique<RooArgSet>();
  RooArgSet(model).snapshot(*_pdfSet, true);
  _pdf = static_cast<RooAbsPdf*>(_pdfSet->find(model.GetName())) ;
 
  _nComp = model._pdfList.size() ;
  _coefThresh.resize(_nComp+1);
  _vars = std::make_unique<RooArgSet>();
  vars.snapshot(*_vars, false);
 
  for (const auto obj : model._pdfList) {
    auto pdf = static_cast<RooAbsPdf*>(obj);
    _gcList.emplace_back(pdf->genContext(vars,prototype,auxProto,verbose));
  }
 
  (static_cast<RooAddModel*>(_pdf))->getProjCache(_vars.get()) ;
  _pdf->recursiveRedirectServers(_theEvent) ;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Attach given set of variables to internal p.d.f. clone
 
void RooAddGenContext::attach(const RooArgSet& args)
{
  _pdf->recursiveRedirectServers(args) ;
 
  // Forward initGenerator call to all components
  for(auto& gc : _gcList) {
    gc->attach(args) ;
  }
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// One-time initialization of generator context. Attach theEvent
/// to internal p.d.f clone and forward initialization call to
/// the component generators
 
void RooAddGenContext::initGenerator(const RooArgSet &theEvent)
{
  _pdf->recursiveRedirectServers(theEvent) ;
 
  if (_isModel) {
    RooAddModel* amod = static_cast<RooAddModel*>(_pdf) ;
    _pcache = amod->getProjCache(_vars.get()) ;
  } else {
    RooAddPdf* apdf = static_cast<RooAddPdf*>(_pdf) ;
    _pcache = apdf->getProjCache(_vars.get(),nullptr) ;
  }
 
  // Forward initGenerator call to all components
  for(auto& gc : _gcList) {
    gc->initGenerator(theEvent) ;
  }
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Randomly choose one of the component contexts to generate this event,
/// with a probability proportional to its coefficient
 
void RooAddGenContext::generateEvent(RooArgSet &theEvent, Int_t remaining)
{
  // Throw a random number to determine which component to generate
  updateThresholds() ;
  double rand = RooRandom::uniform() ;
  for (Int_t i=0 ; i<_nComp ; i++) {
    if (rand>_coefThresh[i] && rand<_coefThresh[i+1]) {
      _gcList[i]->generateEvent(theEvent,remaining) ;
      return ;
    }
  }
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Update the cumulative threshold table from the current coefficient
/// values
 
void RooAddGenContext::updateThresholds()
{
  // Templated lambda to support RooAddModel and RooAddPdf
  auto updateThresholdsImpl = [&](auto* pdf, auto * cache) {
    pdf->updateCoefficients(*cache,_vars.get()) ;
 
    _coefThresh[0] = 0. ;
    for (Int_t i=0 ; i<_nComp ; i++) {
      double coef = pdf->_coefCache[i];
      if(coef < 0.0) {
        std::stringstream errMsgStream;
        errMsgStream << "RooAddGenContext::updateThresholds(): coefficient number " << i << " of the "
                     << pdf->ClassName() << " \"" << pdf->GetName() <<  "\"" << " is negative!"
                     << " The current RooAddGenContext doesn't support negative coefficients."
                     << " Please recreate a new generator context with " << pdf->ClassName() << "::genContext()";
        auto const errMsg = errMsgStream.str();
        cxcoutE(Generation) << errMsg << std::endl;
        throw std::runtime_error(errMsg);
      }
      _coefThresh[i+1] = coef + _coefThresh[i];
    }
  };
 
  _isModel ? updateThresholdsImpl(static_cast<RooAddModel*>(_pdf), _pcache)
           : updateThresholdsImpl(static_cast<RooAddPdf*>(_pdf), _pcache);
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Forward the setProtoDataOrder call to the component generator contexts
 
void RooAddGenContext::setProtoDataOrder(Int_t* lut)
{
  RooAbsGenContext::setProtoDataOrder(lut) ;
  for(auto& gc : _gcList) {
    gc->setProtoDataOrder(lut) ;
  }
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Print the details of the context
 
void RooAddGenContext::printMultiline(std::ostream &os, Int_t content, bool verbose, TString indent) const
{
  RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
  os << indent << "--- RooAddGenContext ---" << std::endl;
  os << indent << "Using PDF ";
  _pdf->printStream(os,kName|kArgs|kClassName,kSingleLine,indent);
 
  os << indent << "List of component generators" << std::endl;
  TString indent2(indent) ;
  indent2.Append("    ") ;
  for(auto& gc : _gcList) {
    gc->printMultiline(os,content,verbose,indent2) ;
  }
}