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 std;

ClassImp(RooStats::HLFactory); ;


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),
    fComboCat(0),
    fComboBkgPdf(0),
    fComboSigBkgPdf(0),
    fComboDataset(0),
    fCombinationDone(false),
    fVerbose(isVerbose),
    fInclusionLevel(0),
    fOwnWs(true){
    TString wsName(name);
    wsName+="_ws";
    fWs = new RooWorkspace(wsName,true);

    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),
    fComboCat(0),
    fComboBkgPdf(0),
    fComboSigBkgPdf(0),
    fComboDataset(0),
    fCombinationDone(false),
    fVerbose(isVerbose),
    fInclusionLevel(0),
    fOwnWs(false){
    fWs=externalWs;
    fSigBkgPdfNames.SetOwner();
    fBkgPdfNames.SetOwner();
    fDatasetsNames.SetOwner();

}

////////////////////////////////////////////////////////////////////////////////

HLFactory::HLFactory():
    TNamed("hlfactory","hlfactory"),
    fComboCat(0),
    fComboBkgPdf(0),
    fComboSigBkgPdf(0),
    fComboDataset(0),
    fCombinationDone(false),
    fVerbose(false),
    fInclusionLevel(0),
    fOwnWs(true){
    fWs = new RooWorkspace("hlfactory_ws",true);

    fSigBkgPdfNames.SetOwner();
    fBkgPdfNames.SetOwner();
    fDatasetsNames.SetOwner();

    }

////////////////////////////////////////////////////////////////////////////////
/// destructor

HLFactory::~HLFactory(){
    if (fComboSigBkgPdf!=NULL)
        delete fComboSigBkgPdf;
    if (fComboBkgPdf!=NULL)
        delete fComboBkgPdf;
    if (fComboDataset!=NULL)
        delete fComboDataset;
    if (fComboCat!=NULL)
        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!=0){
        if (fWs->pdf(SigBkgPdfName)==NULL){
            std::cerr << "Pdf " << SigBkgPdfName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(SigBkgPdfName);
        fSigBkgPdfNames.Add(name);
        }

    if (BkgPdfName!=0){
        if (fWs->pdf(BkgPdfName)==NULL){
            std::cerr << "Pdf " << BkgPdfName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(BkgPdfName);
        fBkgPdfNames.Add(name);
        }

    if (DatasetName!=0){
        if (fWs->data(DatasetName)==NULL){
            std::cerr << "Dataset " << DatasetName << " not found in workspace!\n";
            return -1;
            }
        TObjString* name = new TObjString(DatasetName);
        fDatasetsNames.Add(name);
        }

    if (label!=0){
        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 0;

    if (fComboSigBkgPdf!=NULL)
        return fComboSigBkgPdf;

    if (!fNamesListsConsistent())
        return NULL;

    if (fSigBkgPdfNames.GetSize()==1){
        TString name(((TObjString*)fSigBkgPdfNames.At(0))->String());
        fComboSigBkgPdf=fWs->pdf(name);
        return fComboSigBkgPdf;
        }

    if (!fCombinationDone)
        fCreateCategory();

    RooArgList pdfs("pdfs");

    TIterator* it=fSigBkgPdfNames.MakeIterator();
    TObjString* ostring;
    TObject* obj;
    while ((obj = it->Next())){
        ostring=(TObjString*) obj;
        pdfs.add( *(fWs->pdf(ostring->String())) );
        }
    delete it;

    TString name(GetName());
    name+="_sigbkg";

    TString title(GetName());
    title+="_sigbkg";

    fComboSigBkgPdf=
      new RooSimultaneous(name,
                          title,
                          pdfs,
                          *fComboCat);

    return fComboSigBkgPdf;

}

////////////////////////////////////////////////////////////////////////////////
/// Return the combination of the background only channels.
/// If no background channel is specified a NULL pointer is returned.
/// The factory owns the object.

RooAbsPdf* HLFactory::GetTotBkgPdf(){
    if (fBkgPdfNames.GetSize()==0)
        return 0;

    if (fComboBkgPdf!=NULL)
        return fComboBkgPdf;

    if (!fNamesListsConsistent())
        return NULL;

    if (fBkgPdfNames.GetSize()==1){
        fComboBkgPdf=fWs->pdf(((TObjString*)fBkgPdfNames.First())->String());
        return fComboBkgPdf;
        }

    if (!fCombinationDone)
        fCreateCategory();

    RooArgList pdfs("pdfs");

    TIterator* it = fBkgPdfNames.MakeIterator();
    TObjString* ostring;
    TObject* obj;
    while ((obj = it->Next())){
        ostring=(TObjString*) obj;
        pdfs.add( *(fWs->pdf(ostring->String())) );
        }

    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 NULL pointer is returned.
/// The factory owns the object.

RooDataSet* HLFactory::GetTotDataSet(){
    if (fDatasetsNames.GetSize()==0)
        return 0;

    if (fComboDataset!=NULL)
        return fComboDataset;

    if (!fNamesListsConsistent())
        return NULL;

    if (fDatasetsNames.GetSize()==1){
        fComboDataset=(RooDataSet*)fWs->data(((TObjString*)fDatasetsNames.First())->String());
        return fComboDataset;
        }

    if (!fCombinationDone)
        fCreateCategory();


    TIterator* it = fDatasetsNames.MakeIterator();
    TObjString* ostring;
    TObject* obj = it->Next();
    ostring = (TObjString*) obj;
    fComboDataset = (RooDataSet*) fWs->data(ostring->String()) ;
    if (!fComboDataset) return NULL;
    fComboDataset->Print();
    TString dataname(GetName());
    fComboDataset = new RooDataSet(*fComboDataset,dataname+"_TotData");
    int catindex=0;
    fComboCat->setIndex(catindex);
    fComboDataset->addColumn(*fComboCat);
    while ((obj = it->Next())){
        ostring=(TObjString*) obj;
        catindex++;
        RooDataSet * data = (RooDataSet*)fWs->data(ostring->String());
        if (!data) return NULL;
        RooDataSet* dummy = new RooDataSet(*data,"");
        fComboCat->setIndex(catindex);
        fComboCat->Print();
        dummy->addColumn(*fComboCat);
        fComboDataset->append(*dummy);
        delete dummy;
    }

    delete it;
    return fComboDataset;

}

////////////////////////////////////////////////////////////////////////////////
/// Return the category.
/// The factory owns the object.

RooCategory* HLFactory::GetTotCategory(){
    if (fComboCat!=NULL)
        return fComboCat;

    if (!fNamesListsConsistent())
        return NULL;

    if (!fCombinationDone)
        fCreateCategory();

    return fComboCat;

    }

////////////////////////////////////////////////////////////////////////////////
/// Process an additional configuration file

int HLFactory::ProcessCard(const char* filename){
    return fReadFile(filename,0);
}

////////////////////////////////////////////////////////////////////////////////
/// 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("");

    TObjArray* lines_array = ifileContent.Tokenize("\n");
    TIterator* lineIt=lines_array->MakeIterator();

    bool in_comment=false;
    TString line;
    TObject* line_o;

    while((line_o=(*lineIt)())){ // Start iteration on lines array
        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";
        }

    delete lines_array;
    delete lineIt;

    // Now proceed with the parsing of the stripped file

    lines_array = ifileContentStripped.Tokenize(";");
    lineIt=lines_array->MakeIterator();
    in_comment=false;

    const int nNeutrals=2;
    TString neutrals[nNeutrals]={"\t"," "};

    while((line_o=(*lineIt)())){

        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);
        }

    delete lineIt;
    delete lines_array;

    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);

    TIterator* it=fLabelsNames.MakeIterator();
    TObjString* ostring;
    TObject* obj;
    while ((obj = it->Next())){
        ostring=(TObjString*) obj;
        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);
      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());

        TFile* ifile=TFile::Open(rootfile_name);
        if (ifile==0)
            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);
          }
        delete ifile;
        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);

        return 0;
        }

    else { // In case we do not know what to do we pipe it..
        fWs->factory(line);
        }

    return 0;

}