HLFactory.cxx

Go to the documentation of this file.

// @(#)root/roostats:$Id$
// Author: Danilo Piparo   25/08/2009
 
 
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
 
////////////////////////////////////////////////////////////////////////////////
 
 
#include <iostream>
#include <fstream>
 
#include "RooStats/HLFactory.h"
#include "TFile.h"
#include "TObject.h"
#include "TObjArray.h"
#include "TObjString.h"
 
#include "RooSimultaneous.h"
 
/** \class RooStats::HLFactory
    \ingroup Roostats
 
HLFactory is an High Level model Factory allows you to
describe your models in a configuration file
(_datacards_) acting as an interface with the RooFactoryWSTool.
Moreover it provides tools for the combination of models and datasets.
 
*/
 
 
 
using namespace RooStats;
using namespace RooFit;
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor with the name of the config file to interpret and the
/// verbosity flag. The extension for the config files is assumed to
/// be ".rs".
 
HLFactory::HLFactory(const char *name, const char *fileName, bool isVerbose)
   : TNamed(name, name), fVerbose(isVerbose), fOwnWs(true)
{
   TString wsName(name);
   wsName += "_ws";
   fWs = new RooWorkspace(wsName);
 
   fSigBkgPdfNames.SetOwner();
   fBkgPdfNames.SetOwner();
   fDatasetsNames.SetOwner();
 
   // Start the parsing
   fReadFile(fileName);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor without a card but with an external workspace.
 
HLFactory::HLFactory(const char *name, RooWorkspace *externalWs, bool isVerbose)
   : TNamed(name, name), fVerbose(isVerbose), fWs(externalWs)
{
   fSigBkgPdfNames.SetOwner();
   fBkgPdfNames.SetOwner();
   fDatasetsNames.SetOwner();
}
 
////////////////////////////////////////////////////////////////////////////////
 
HLFactory::HLFactory() : TNamed("hlfactory", "hlfactory"), fWs(new RooWorkspace("hlfactory_ws")), fOwnWs(true)
{
   fSigBkgPdfNames.SetOwner();
   fBkgPdfNames.SetOwner();
   fDatasetsNames.SetOwner();
}
 
////////////////////////////////////////////////////////////////////////////////
/// destructor
 
HLFactory::~HLFactory(){
    if (fComboSigBkgPdf!=nullptr)
        delete fComboSigBkgPdf;
    if (fComboBkgPdf!=nullptr)
        delete fComboBkgPdf;
    if (fComboDataset!=nullptr)
        delete fComboDataset;
    if (fComboCat!=nullptr)
        delete fComboCat;
 
    if (fOwnWs)
        delete fWs;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Add a channel to the combination. The channel can be specified as:
///  - A signal plus background pdf
///  - A background only pdf
///  - A dataset
/// Once the combination of the pdfs is done, no more channels should be
/// added.
 
int HLFactory::AddChannel(const char* label,
                          const char* SigBkgPdfName,
                          const char* BkgPdfName,
                          const char* DatasetName){
    if (fCombinationDone){
        std::cerr << "Cannot add anymore channels. "
                  << "Combination already carried out.\n";
        return -1;
        }
 
    if (SigBkgPdfName!=nullptr){
        if (fWs->pdf(SigBkgPdfName)==nullptr){
            std::cerr << "Pdf " << SigBkgPdfName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(SigBkgPdfName);
        fSigBkgPdfNames.Add(name);
        }
 
    if (BkgPdfName!=nullptr){
        if (fWs->pdf(BkgPdfName)==nullptr){
            std::cerr << "Pdf " << BkgPdfName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(BkgPdfName);
        fBkgPdfNames.Add(name);
        }
 
    if (DatasetName!=nullptr){
        if (fWs->data(DatasetName)==nullptr){
            std::cerr << "Dataset " << DatasetName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(DatasetName);
        fDatasetsNames.Add(name);
        }
 
    if (label!=nullptr){
        TObjString* name = new TObjString(label);
        fLabelsNames.Add(name);
        }
    return 0;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the combination of the signal plus background channels.
/// The factory owns the object.
 
RooAbsPdf* HLFactory::GetTotSigBkgPdf(){
    if (fSigBkgPdfNames.GetSize()==0)
        return nullptr;
 
    if (fComboSigBkgPdf!=nullptr)
        return fComboSigBkgPdf;
 
    if (!fNamesListsConsistent())
        return nullptr;
 
    if (fSigBkgPdfNames.GetSize()==1){
        fComboSigBkgPdf=fWs->pdf(static_cast<TObjString*>(fSigBkgPdfNames.At(0))->String().Data());
        return fComboSigBkgPdf;
        }
 
    if (!fCombinationDone)
        fCreateCategory();
 
    RooArgList pdfs("pdfs");
 
    for(auto * ostring : static_range_cast<TObjString*>(fSigBkgPdfNames)) {
        pdfs.add( *(fWs->pdf(ostring->String().Data())) );
    }
 
    std::string name(GetName());
    name+="_sigbkg";
 
    std::string title(GetName());
    title+="_sigbkg";
 
    fComboSigBkgPdf=
      new RooSimultaneous(name.c_str(),
                          title.c_str(),
                          pdfs,
                          *fComboCat);
 
    return fComboSigBkgPdf;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the combination of the background only channels.
/// If no background channel is specified a nullptr pointer is returned.
/// The factory owns the object.
 
RooAbsPdf* HLFactory::GetTotBkgPdf(){
    if (fBkgPdfNames.GetSize()==0)
        return nullptr;
 
    if (fComboBkgPdf!=nullptr)
        return fComboBkgPdf;
 
    if (!fNamesListsConsistent())
        return nullptr;
 
    if (fBkgPdfNames.GetSize()==1){
        fComboBkgPdf=fWs->pdf(static_cast<TObjString*>(fBkgPdfNames.First())->String().Data());
        return fComboBkgPdf;
        }
 
    if (!fCombinationDone)
        fCreateCategory();
 
    RooArgList pdfs("pdfs");
 
    for(auto * ostring : static_range_cast<TObjString*>(fBkgPdfNames)) {
        pdfs.add( *fWs->pdf(ostring->String().Data()) );
    }
 
    TString name(GetName());
    name+="_bkg";
 
    TString title(GetName());
    title+="_bkg";
 
    fComboBkgPdf=
      new RooSimultaneous(name,
                          title,
                          pdfs,
                          *fComboCat);
 
    return fComboBkgPdf;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the combination of the datasets.
/// If no dataset is specified a nullptr pointer is returned.
/// The factory owns the object.
 
RooDataSet* HLFactory::GetTotDataSet(){
    if (fDatasetsNames.GetSize()==0)
        return nullptr;
 
    if (fComboDataset!=nullptr)
        return fComboDataset;
 
    if (!fNamesListsConsistent())
        return nullptr;
 
    if (fDatasetsNames.GetSize()==1){
        fComboDataset=static_cast<RooDataSet*>(fWs->data(static_cast<TObjString*>(fDatasetsNames.First())->String().Data()));
        return fComboDataset;
        }
 
    if (!fCombinationDone)
        fCreateCategory();
 
 
    auto it = fDatasetsNames.begin();
    TObjString* ostring;
    ostring = static_cast<TObjString*>(*it);
    ++it;
    fComboDataset = static_cast<RooDataSet*>(fWs->data(ostring->String().Data())) ;
    if (!fComboDataset) return nullptr;
    fComboDataset->Print();
    TString dataname(GetName());
    fComboDataset = new RooDataSet(*fComboDataset,dataname+"_TotData");
    int catindex=0;
    fComboCat->setIndex(catindex);
    fComboDataset->addColumn(*fComboCat);
    for(; it != fDatasetsNames.end() ; ++it) {
        ostring = static_cast<TObjString*>(*it);
        catindex++;
        RooDataSet * data = static_cast<RooDataSet*>(fWs->data(ostring->String().Data()));
        if (!data) return nullptr;
        RooDataSet* dummy = new RooDataSet(*data,"");
        fComboCat->setIndex(catindex);
        fComboCat->Print();
        dummy->addColumn(*fComboCat);
        fComboDataset->append(*dummy);
        delete dummy;
    }
 
    return fComboDataset;
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the category.
/// The factory owns the object.
 
RooCategory* HLFactory::GetTotCategory(){
    if (fComboCat!=nullptr)
        return fComboCat;
 
    if (!fNamesListsConsistent())
        return nullptr;
 
    if (!fCombinationDone)
        fCreateCategory();
 
    return fComboCat;
 
    }
 
////////////////////////////////////////////////////////////////////////////////
/// Process an additional configuration file
 
int HLFactory::ProcessCard(const char* filename){
    return fReadFile(filename,false);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Parses the configuration file. The objects can be specified following
/// the rules of the RooFactoryWSTool, plus some more flexibility.
///
/// The official format for the datacards is ".rs".
///
/// All the instructions end with a ";" (like in C++).
///
/// Carriage returns and white lines are irrelevant but advised since they
/// improve readability (like in C++).
///
/// The `(Roo)ClassName::objname(description)` can be replaced with the more
/// "pythonic" `objname = (Roo)ClassName(description)`.
///
/// The comments can be specified with a "//" if on a single line or with
/// "multiple lines" in C/C++ like comments.
///
/// The `"#include path/to/file.rs"` statement triggers the inclusion of a
/// configuration fragment.
///
/// The `"import myobject:myworkspace:myrootfile"` will add to the Workspace
/// the object myobject located in myworkspace recorded in myrootfile.
/// Alternatively, one could choose the `"import myobject:myrootfile"` in case
/// no Workspace is present.
///
/// The `"echo"` statement prompts a message on screen.
 
int HLFactory::fReadFile(const char*fileName, bool is_included){
    // Check the deepness of the inclusion
    if (is_included) {
        fInclusionLevel+=1;
    } else {
        fInclusionLevel = 0;
    }
 
    const int maxDeepness=50;
    if (fInclusionLevel>maxDeepness){
        TString warning("The inclusion stack is deeper than ");
        warning+=maxDeepness;
        warning+=". Is this a recursive inclusion?";
        Warning("fReadFile", "%s", warning.Data());
        }
 
 
    // open the config file and go through it
    std::ifstream ifile(fileName);
 
    if(ifile.fail()){
        TString error("File ");
        error+=fileName;
        error+=" could not be opened.";
        Error("fReadFile", "%s", error.Data());
        return -1;
        }
 
    TString ifileContent("");
    ifileContent.ReadFile(ifile);
    ifile.close();
 
    // Tokenise the file using the "\n" char and parse it line by line to strip
    // the comments.
    TString ifileContentStripped("");
 
    std::unique_ptr<TObjArray> lines_array{ifileContent.Tokenize("\n")};
 
    bool in_comment=false;
 
    // Start iteration on lines array
    for(TObject * line_o : *lines_array) {
        TString line = (static_cast<TObjString*>(line_o))->GetString();
 
        // Are we in a multiline comment?
        if (in_comment) {
            if (line.EndsWith("*/")){
                in_comment=false;
                if (fVerbose) Info("fReadFile","Out of multiline comment ...");
 
                continue;
            }
        }
 
        // Was line a single line comment?
 
        if ((line.BeginsWith("/*") && line.EndsWith("*/")) ||
            line.BeginsWith("//")){
            if (fVerbose) Info("fReadFile","In single line comment ...");
            continue;
            }
 
        // Did a multiline comment just begin?
        if (line.BeginsWith("/*")){
            in_comment=true;
            if (fVerbose) Info("fReadFile","In multiline comment ...");
            continue;
            }
 
        ifileContentStripped+=line+"\n";
        }
 
    // Now proceed with the parsing of the stripped file
 
    lines_array.reset(ifileContentStripped.Tokenize(";"));
    in_comment=false;
 
    const int nNeutrals=2;
    TString neutrals[nNeutrals]={"\t"," "};
 
    for(TObject * line_o : *lines_array) {
 
        TString line = (static_cast<TObjString*>(line_o))->GetString();
 
        // Strip spaces at the beginning and the end of the line
        line.Strip(TString::kBoth,' ');
 
        // Put the single statement in one single line
        line.ReplaceAll("\n","");
 
        // Do we have an echo statement? "A la RooFit"
        if (line.BeginsWith("echo")){
            line = line(5,line.Length()-1);
            if (fVerbose)
              std::cout << "Echoing line " << line.Data() << std::endl;
            std::cout << "[" << GetName() << "] echo: "
                      << line.Data() << std::endl;
            continue;
            }
 
        // Spaces and tabs at this point are not needed.
        for (int i=0;i<nNeutrals;++i)
            line.ReplaceAll(neutrals[i],"");
 
 
        if (fVerbose) Info("fReadFile","Reading --> %s <--", line.Data());
 
        // Was line a white space?
        if (line == ""){
            if (fVerbose) Info("fReadFile", "%s", "Empty line: skipping ...");
            continue;
            }
 
        // Do we have an include statement?
        // We treat this recursively.
        if (line.BeginsWith("#include")){
            line.ReplaceAll("#include","");
            if (fVerbose) Info("fReadFile","Reading included file...");
            fReadFile(line,true);
            continue;
            }
 
        // We parse the line
        if (fVerbose) Info("fReadFile","Parsing the line...");
        fParseLine(line);
        }
 
    return 0;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Builds the category necessary for the mutidimensional models. Its name
/// will be `<HLFactory name>_category` and the types are specified by the
/// model labels.
 
void HLFactory::fCreateCategory(){
    fCombinationDone=true;
 
    TString name(GetName());
    name+="_category";
 
    TString title(GetName());
    title+="_category";
 
    fComboCat=new RooCategory(name,title);
 
    for (auto * ostring : static_range_cast<TObjString*>(fLabelsNames)) {
        fComboCat->defineType(ostring->String());
        }
 
    }
 
////////////////////////////////////////////////////////////////////////////////
/// Check the number of entries in each list. If not the same and the list
/// is not empty prompt an error.
 
bool HLFactory::fNamesListsConsistent(){
        if ((fSigBkgPdfNames.GetEntries() == fBkgPdfNames.GetEntries() || fBkgPdfNames.GetEntries() == 0) &&
            (fSigBkgPdfNames.GetEntries() == fDatasetsNames.GetEntries() || fDatasetsNames.GetEntries() == 0) &&
            (fSigBkgPdfNames.GetEntries() == fLabelsNames.GetEntries() || fLabelsNames.GetEntries() == 0)) {
        return true;
        } else {
        std::cerr << "The number of datasets and models added as channels "
                  << " is not the same!\n";
        return false;
        }
    }
 
////////////////////////////////////////////////////////////////////////////////
/// Parse a single line and puts the content in the RooWorkSpace
 
int HLFactory::fParseLine(TString& line){
    if (fVerbose) Info("fParseLine", "Parsing line: %s", line.Data());
 
    TString new_line("");
 
    const int nequals = line.CountChar('=');
 
    // Build with the factory a var or cat, or pipe the command directly.
 
    if (line.Contains("::") || // It is a ordinary statement
        nequals==0 || //it is a RooRealVar or cat with 0,1,2,3.. indexes
        (line.Contains("[") &&
         line.Contains("]") &&
         nequals>0 &&    // It is a cat like "tag[B0=1,B0bar=-1]"
         ! line.Contains("(") &&
         ! line.Contains(")"))) {
      fWs->factory(line.Data());
      return 0;
      }
 
    // Transform the line o_name = o_class(o_descr) in o_class::o_name(o_descr)
    if (nequals==1 ||
        (nequals > 1 &&  line.Contains("SIMUL"))){
 
        // Divide the line in 3 components: o_name,o_class and o_descr
        // assuming that o_name=o_class(o_descr)
        const int equal_index=line.First('=');
        const int par_index=line.First('(');
        TString o_name(line(0,equal_index));
        TString o_class(line(equal_index+1,par_index-equal_index-1));
        TString o_descr(line(par_index+1,line.Length()-par_index-2));
 
        if (fVerbose) Info("fParseLine", "o_name=%s o_class=%s o_descr=%s",
                           o_name.Data(), o_class.Data(), o_descr.Data());
 
        // Now two cases either we wanna produce an object or import something
        // under a new name.
        if (o_class =="import"){// import a generic TObject into the WS
        // Now see if we have a workspace or not, according to the number of
        // entries in the description..
 
        TObjArray* descr_array = o_descr.Tokenize(",");
 
        const int n_descr_parts=descr_array->GetEntries();
 
        if (n_descr_parts<2 || n_descr_parts>3)
          Error("fParseLine","Import wrong syntax: cannot process %s", o_descr.Data());
 
        TString obj_name (static_cast<TObjString*>(descr_array->At(n_descr_parts-1))->GetString());
        TString ws_name("");
        TString rootfile_name (static_cast<TObjString*>(descr_array->At(0))->GetString());
 
        std::unique_ptr<TFile> ifile{TFile::Open(rootfile_name)};
        if (ifile==nullptr)
            return 1;
 
        if (n_descr_parts==3){// in presence of a Ws
          o_descr.ReplaceAll(",",":");
          fWs->import(o_descr);
          }
        else if(n_descr_parts==2){ // in presence of an object in rootfile
          if (fVerbose) {
             Info("fParseLine", "Importing %s from %s under the name of %s", obj_name.Data(), rootfile_name.Data(),
                  o_name.Data());
          }
          TObject* the_obj=ifile->Get(obj_name);
          fWs->import(*the_obj,o_name);
          }
        return 0;
        } // end of import block
 
        new_line=o_class+"::"+o_name+"("+o_descr+")";
 
        if (fVerbose){
            std::cout << "DEBUG: line: " << line.Data() << std::endl;
            std::cout << "DEBUG: new_line: " << new_line.Data() << std::endl;
            }
 
        fWs->factory(new_line.Data());
 
        return 0;
        }
 
    else { // In case we do not know what to do we pipe it..
        fWs->factory(line.Data());
        }
 
    return 0;
 
}