Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
171 changes: 83 additions & 88 deletions roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@
#include <fstream>
#include <iomanip>
#include <memory>
#include <set>
#include <utility>

constexpr double alphaLow = -5.0;
Expand All @@ -94,17 +95,17 @@ std::vector<double> histToVector(TH1 const &hist)

// use this order for safety on library loading
using namespace RooStats;
using std::string, std::vector, std::make_unique, std::pair, std::unique_ptr, std::map;
using std::string, std::vector;

using namespace RooStats::HistFactory::Detail;
using namespace RooStats::HistFactory::Detail::MagicConstants;

namespace RooStats::HistFactory {

namespace RooStats{
namespace HistFactory{

HistoToWorkspaceFactoryFast::HistoToWorkspaceFactoryFast(RooStats::HistFactory::Measurement& measurement) :
HistoToWorkspaceFactoryFast{measurement, Configuration{}} {}
HistoToWorkspaceFactoryFast::HistoToWorkspaceFactoryFast(RooStats::HistFactory::Measurement &measurement)
: HistoToWorkspaceFactoryFast{measurement, Configuration{}}
{
}

HistoToWorkspaceFactoryFast::HistoToWorkspaceFactoryFast(RooStats::HistFactory::Measurement& measurement,
Configuration const& cfg) :
Expand All @@ -130,9 +131,8 @@ namespace HistFactory{
// Make a ModelConfig and configure it
ModelConfig * proto_config = static_cast<ModelConfig *>(ws_single->obj("ModelConfig"));
if( proto_config == nullptr ) {
std::cout << "Error: Did not find 'ModelConfig' object in file: " << ws_single->GetName()
<< std::endl;
throw hf_exc();
cxcoutFHF << "Error: Did not find 'ModelConfig' object in file: " << ws_single->GetName() << std::endl;
throw hf_exc();
}

if( measurement.GetPOIList().empty() ) {
Expand All @@ -153,9 +153,9 @@ namespace HistFactory{
params.add(*poi);
}
else {
std::cout << "WARNING: Can't find parameter of interest: " << poi_name
<< " in Workspace. Not setting in ModelConfig." << std::endl;
//throw hf_exc();
cxcoutWHF << "WARNING: Can't find parameter of interest: " << poi_name
<< " in Workspace. Not setting in ModelConfig." << std::endl;
// throw hf_exc();
}
}
proto_config->SetParametersOfInterest(params);
Expand Down Expand Up @@ -196,9 +196,8 @@ namespace HistFactory{
cxcoutPHF << "Importing Asimov dataset" << std::endl;
bool failure = ws_single->import(*asimov_dataset, RooFit::Rename(AsimovName.c_str()));
if( failure ) {
std::cout << "Error: Failed to import Asimov dataset: " << AsimovName
<< std::endl;
throw hf_exc();
cxcoutFHF << "Error: Failed to import Asimov dataset: " << AsimovName << std::endl;
throw hf_exc();
}

// Load the snapshot at the end of every loop iteration
Expand Down Expand Up @@ -229,9 +228,9 @@ namespace HistFactory{
// Create a workspace for a SingleChannel from the Measurement Object
std::unique_ptr<RooWorkspace> ws_single{this->MakeSingleChannelWorkspace(measurement, channel)};
if( ws_single == nullptr ) {
cxcoutF(HistFactory) << "Error: Failed to make Single-Channel workspace for channel: " << ch_name
<< " and measurement: " << measurement.GetName() << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Failed to make Single-Channel workspace for channel: " << ch_name
<< " and measurement: " << measurement.GetName() << std::endl;
throw hf_exc();
}

// Finally, configure that workspace based on
Expand Down Expand Up @@ -277,8 +276,8 @@ namespace HistFactory{
HistoToWorkspaceFactoryFast histFactory(measurement, config);

// Loop over the channels and create the individual workspaces
vector<std::unique_ptr<RooWorkspace>> channel_workspaces;
vector<string> channel_names;
std::vector<std::unique_ptr<RooWorkspace>> channel_workspaces;
std::vector<std::string> channel_names;

for(HistFactory::Channel& channel : measurement.GetChannels()) {

Expand Down Expand Up @@ -748,9 +747,9 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
}

if( ! channel.CheckHistograms() ) {
std::cout << "MakeSingleChannelWorkspace: Channel: " << channel.GetName()
<< " has uninitialized histogram pointers" << std::endl;
throw hf_exc();
cxcoutFHF << "MakeSingleChannelWorkspace: Channel: " << channel.GetName()
<< " has uninitialized histogram pointers" << std::endl;
throw hf_exc();
}


Expand Down Expand Up @@ -796,7 +795,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
//
auto protoOwner = std::make_unique<RooWorkspace>(channel_name.c_str(), (channel_name+" workspace").c_str());
RooWorkspace &proto = *protoOwner;
auto proto_config = make_unique<ModelConfig>("ModelConfig", &proto);
auto proto_config = std::make_unique<ModelConfig>("ModelConfig", &proto);

// preprocess functions
for(auto const& func : fPreprocessFunctions){
Expand All @@ -814,8 +813,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
std::vector<std::vector<RooAbsArg*>> allSampleHistFuncs;
std::vector<RooProduct*> sampleScaleFactors;

std::vector< pair<string,string> > statNamePairs;
std::vector< pair<const TH1*, std::unique_ptr<TH1>> > statHistPairs; // <nominal, error>
std::vector<std::pair<string, string>> statNamePairs;
std::vector<std::pair<const TH1 *, std::unique_ptr<TH1>>> statHistPairs; // <nominal, error>
const std::string statFuncName = "mc_stat_" + channel_name;

string prefix;
Expand Down Expand Up @@ -918,9 +917,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
if( sample.GetStatError().GetActivate() ) {

if( fObsNameVec.size() > 3 ) {
cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
<< std::endl;
throw hf_exc();
cxcoutFHF << "Cannot include Stat Error for histograms of more than 3 dimensions." << std::endl;
throw hf_exc();
} else {

// If we are using StatUncertainties, we multiply this object
Expand Down Expand Up @@ -1018,9 +1016,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
if( !sample.GetShapeFactorList().empty() ) {

if( fObsNameVec.size() > 3 ) {
cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
<< std::endl;
throw hf_exc();
cxcoutFHF << "Cannot include Stat Error for histograms of more than 3 dimensions." << std::endl;
throw hf_exc();
} else {

cxcoutI(HistFactory) << "Sample: " << sample.GetName() << " in channel: " << channel_name
Expand Down Expand Up @@ -1086,9 +1083,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
if( !sample.GetShapeSysList().empty() ) {

if( fObsNameVec.size() > 3 ) {
cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
<< std::endl;
throw hf_exc();
cxcoutFHF << "Cannot include Stat Error for histograms of more than 3 dimensions.\n";
throw hf_exc();
}

// List of ShapeSys ParamHistFuncs
Expand Down Expand Up @@ -1210,11 +1206,11 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo

// Create the histogram of (binwise)
// stat uncertainties:
unique_ptr<TH1> fracStatError( MakeScaledUncertaintyHist( channel_name + "_StatUncert" + "_RelErr", statHistPairs) );
std::unique_ptr<TH1> fracStatError(
MakeScaledUncertaintyHist(channel_name + "_StatUncert" + "_RelErr", statHistPairs));
if( fracStatError == nullptr ) {
cxcoutE(HistFactory) << "Error: Failed to make ScaledUncertaintyHist for: "
<< channel_name + "_StatUncert" + "_RelErr" << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Failed to make ScaledUncertaintyHist for: " << channel_name + "_StatUncert" + "_RelErr\n";
throw hf_exc();
}

// Using this TH1* of fractinal stat errors,
Expand Down Expand Up @@ -1281,19 +1277,19 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
for(unsigned int i=0; i<constraintTermNames.size(); ++i){
RooAbsArg* proto_arg = proto.arg(constraintTermNames[i]);
if( proto_arg==nullptr ) {
cxcoutF(HistFactory) << "Error: Cannot find arg set: " << constraintTermNames.at(i)
<< " in workspace: " << proto.GetName() << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Cannot find arg set: " << constraintTermNames.at(i)
<< " in workspace: " << proto.GetName() << std::endl;
throw hf_exc();
}
constraintTerms.add( *proto_arg );
// constraintTerms.add(* proto_arg(proto.arg(constraintTermNames[i].c_str())) );
}
for(unsigned int i=0; i<likelihoodTermNames.size(); ++i){
RooAbsArg* proto_arg = (proto.arg(likelihoodTermNames[i]));
if( proto_arg==nullptr ) {
cxcoutF(HistFactory) << "Error: Cannot find arg set: " << likelihoodTermNames.at(i)
<< " in workspace: " << proto.GetName() << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Cannot find arg set: " << likelihoodTermNames.at(i)
<< " in workspace: " << proto.GetName() << std::endl;
throw hf_exc();
}
likelihoodTerms.add( *proto_arg );
}
Expand Down Expand Up @@ -1322,10 +1318,10 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
<< "\timport model into workspace"
<< "\n-----------------------------------------\n" << std::endl;

auto model = make_unique<RooProdPdf>(
("model_"+channel_name).c_str(), // MB : have changed this into conditional pdf. Much faster for toys!
"product of Poissons across bins for a single channel",
constraintTerms, RooFit::Conditional(likelihoodTerms,observables));
auto model = std::make_unique<RooProdPdf>(
("model_" + channel_name).c_str(), // MB : have changed this into conditional pdf. Much faster for toys!
"product of Poissons across bins for a single channel", constraintTerms,
RooFit::Conditional(likelihoodTerms, observables));
// can give channel a title by setting title of corresponding data histogram
if (channel.GetData().GetHisto() && strlen(channel.GetData().GetHisto()->GetTitle())>0) {
model->SetTitle(channel.GetData().GetHisto()->GetTitle());
Expand Down Expand Up @@ -1355,7 +1351,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
// actually create the files with the stored per-channel workspaces.
// Otherwise, we just spend time calculating something that gets thrown
// away anyway (for the combined workspace, we'll create a new Asimov).
unique_ptr<RooAbsData> asimov_dataset(AsymptoticCalculator::GenerateAsimovData(*model, observables));
std::unique_ptr<RooAbsData> asimov_dataset(AsymptoticCalculator::GenerateAsimovData(*model, observables));
proto.import(*asimov_dataset, RooFit::Rename("asimovData"));
}

Expand All @@ -1378,17 +1374,17 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo

for(auto const& data : channel.GetAdditionalData()) {
if(data.GetName().empty()) {
cxcoutF(HistFactory) << "Error: Additional Data histogram for channel: " << channel.GetName()
<< " has no name! The name always needs to be set for additional datasets, "
<< "either via the \"Name\" tag in the XML or via RooStats::HistFactory::Data::SetName()." << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Additional Data histogram for channel: " << channel.GetName()
<< " has no name! The name always needs to be set for additional datasets, "
<< "either via the \"Name\" tag in the XML or via RooStats::HistFactory::Data::SetName().\n";
throw hf_exc();
}
std::string const& dataName = data.GetName();
TH1 const* mnominal = data.GetHisto();
if( !mnominal ) {
cxcoutF(HistFactory) << "Error: Additional Data histogram for channel: " << channel.GetName()
<< " with name: " << dataName << " is nullptr" << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Additional Data histogram for channel: " << channel.GetName()
<< " with name: " << dataName << " is nullptr\n";
throw hf_exc();
}

// THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
Expand Down Expand Up @@ -1483,18 +1479,23 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
Error("MakeCombinedModel","Input vector of workspace has an invalid size - return a nullptr");
return nullptr;
}
std::set<std::string> ch_names_set{ch_names.begin(), ch_names.end()};
if (ch_names.size() != ch_names_set.size()) {
Error("MakeCombinedModel", "Input vector of channel names has duplicate names - return a nullptr");
return nullptr;
}

//
/// These things were used for debugging. Maybe useful in the future
//

map<string, RooAbsPdf*> pdfMap;
vector<RooAbsPdf*> models;
std::map<string, RooAbsPdf *> pdfMap;
vector<RooAbsPdf *> models;

RooArgList obsList;
for(unsigned int i = 0; i< ch_names.size(); ++i){
obsList.add(*static_cast<ModelConfig *>(chs[i]->obj("ModelConfig"))->GetObservables());
}
RooArgList obsList;
for (unsigned int i = 0; i < ch_names.size(); ++i) {
obsList.add(*static_cast<ModelConfig *>(chs[i]->obj("ModelConfig"))->GetObservables());
}
cxcoutI(HistFactory) <<"full list of observables:\n" << obsList << std::endl;

RooArgSet globalObs;
Expand All @@ -1512,8 +1513,10 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
RooWorkspace * ch=chs[i].get();

RooAbsPdf* model = ch->pdf("model_"+channel_name);
if(!model) std::cout <<"failed to find model for channel"<< std::endl;
// std::cout << "int = " << model->createIntegral(*obsN)->getVal() << std::endl;
if (!model) {
cxcoutFHF << "failed to find model for channel\n";
throw hf_exc();
}
models.push_back(model);
globalObs.add(*ch->set("globalObservables"), /*silent=*/true); // silent because observables might exist in other channel.

Expand Down Expand Up @@ -1621,8 +1624,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
combined->import( *asimov_combined, RooFit::Rename("asimovData"));
}
else {
std::cout << "Error: Failed to create combined asimov dataset" << std::endl;
throw hf_exc();
cxcoutFHF << "Error: Failed to create combined asimov dataset\n";
throw hf_exc();
}

return RooFit::makeOwningPtr(std::move(combined));
Expand Down Expand Up @@ -1654,23 +1657,19 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo

// Check that histError != NAN
if( histError != histError ) {
std::cout << "Warning: In histogram " << Nominal->GetName()
<< " bin error for bin " << i_bin
<< " is NAN. Not using Error!!!"
<< std::endl;
throw hf_exc();
//histError = sqrt( histContent );
//histError = 0;
cxcoutFHF << "Warning: In histogram " << Nominal->GetName() << " bin error for bin " << i_bin
<< " is NAN. Not using Error!!!\n";
throw hf_exc();
// histError = sqrt( histContent );
// histError = 0;
}

// Check that histError ! < 0
if( histError < 0 ) {
std::cout << "Warning: In histogram " << Nominal->GetName()
<< " bin error for bin " << binNumber
<< " is < 0. Setting Error to 0"
<< std::endl;
//histError = sqrt( histContent );
histError = 0;
cxcoutWHF << "Warning: In histogram " << Nominal->GetName() << " bin error for bin " << binNumber
<< " is < 0. Setting Error to 0" << std::endl;
// histError = sqrt( histContent );
histError = 0;
}

ErrorHist->SetBinContent( binNumber, histError );
Expand Down Expand Up @@ -1743,10 +1742,9 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
double histError = error->GetBinContent( binNumber );

if( histError != histError ) {
cxcoutE(HistFactory) << "In histogram " << error->GetName()
<< " bin error for bin " << binNumber
<< " is NAN. Not using error!!";
throw hf_exc();
cxcoutFHF << "In histogram " << error->GetName() << " bin error for bin " << binNumber
<< " is NAN. Not using error!!";
throw hf_exc();
}

TotalBinContent.at(i_bins) += histValue;
Expand Down Expand Up @@ -1812,7 +1810,4 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
return std::unique_ptr<TH1>(ErrorHist);
}


} // namespace RooStats
} // namespace HistFactory

} // namespace RooStats::HistFactory
Loading