// @(#)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. * *************************************************************************/ //___________________________________________________ /* BEGIN_HTML

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.

END_HTML */ #include #include #include "RooStats/HLFactory.h" #include "TFile.h" #include "TObject.h" #include "TObjArray.h" #include "TObjString.h" #include "RooSimultaneous.h" using namespace std; ClassImp(RooStats::HLFactory) ; using namespace RooStats; using namespace RooFit; //_______________________________________________________ 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){ // 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". TString wsName(name); wsName+="_ws"; fWs = new RooWorkspace(wsName,true); fSigBkgPdfNames.SetOwner(); fBkgPdfNames.SetOwner(); fDatasetsNames.SetOwner(); // Start the parsing fReadFile(fileName); } //_______________________________________________________ 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){ // Constructor without a card but with an exrernal workspace. 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(); } //_______________________________________________________ HLFactory::~HLFactory(){ // destructor if (fComboSigBkgPdf!=NULL) delete fComboSigBkgPdf; if (fComboBkgPdf!=NULL) delete fComboBkgPdf; if (fComboDataset!=NULL) delete fComboDataset; if (fComboCat!=NULL) delete fComboCat; if (fOwnWs) delete fWs; } //_______________________________________________________ int HLFactory::AddChannel(const char* label, const char* SigBkgPdfName, const char* BkgPdfName, const char* DatasetName){ // 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. 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; } //_______________________________________________________ RooAbsPdf* HLFactory::GetTotSigBkgPdf(){ // Return the combination of the signal plus background channels. // The facory owns the object. 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; } //_______________________________________________________ RooAbsPdf* HLFactory::GetTotBkgPdf(){ // Return the combination of the background only channels. // If no background channel is specified a NULL pointer is returned. // The facory owns the object. 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; } //_______________________________________________________ RooDataSet* HLFactory::GetTotDataSet(){ // Return the combination of the datasets. // If no dataset is specified a NULL pointer is returned. // The facory owns the object. 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; } //_______________________________________________________ RooCategory* HLFactory::GetTotCategory(){ // Return the category. // The facory owns the object. if (fComboCat!=NULL) return fComboCat; if (!fNamesListsConsistent()) return NULL; if (!fCombinationDone) fCreateCategory(); return fComboCat; } //_______________________________________________________ int HLFactory::ProcessCard(const char* filename){ // Process an additional configuration file return fReadFile(filename,0); } //_______________________________________________________ int HLFactory::fReadFile(const char*fileName, bool is_included){ // 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 adviced 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 // /* */ if on multiple lines (like in C++). // // 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. // 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(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(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 %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; } //_______________________________________________________ void HLFactory::fCreateCategory(){ // Builds the category necessary for the mutidimensional models. Its name // will be _category and the types are specified by the // model labels. 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()); } } //_______________________________________________________ bool HLFactory::fNamesListsConsistent(){ // Check the number of entries in each list. If not the same and the list // is not empty prompt an error. 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; } } //_______________________________________________________ int HLFactory::fParseLine(TString& line){ // Parse a single line and puts the content in the RooWorkSpace 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(descr_array->At(n_descr_parts-1))->GetString()); TString ws_name(""); TString rootfile_name (static_cast(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; }