diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx index f5d9a473ffcf5..7033943d19133 100644 --- a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx +++ b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx @@ -71,6 +71,7 @@ #include #include #include +#include #include constexpr double alphaLow = -5.0; @@ -94,17 +95,17 @@ std::vector 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) : @@ -130,9 +131,8 @@ namespace HistFactory{ // Make a ModelConfig and configure it ModelConfig * proto_config = static_cast(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() ) { @@ -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); @@ -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 @@ -229,9 +228,9 @@ namespace HistFactory{ // Create a workspace for a SingleChannel from the Measurement Object std::unique_ptr 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 @@ -277,8 +276,8 @@ namespace HistFactory{ HistoToWorkspaceFactoryFast histFactory(measurement, config); // Loop over the channels and create the individual workspaces - vector> channel_workspaces; - vector channel_names; + std::vector> channel_workspaces; + std::vector channel_names; for(HistFactory::Channel& channel : measurement.GetChannels()) { @@ -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(); } @@ -796,7 +795,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo // auto protoOwner = std::make_unique(channel_name.c_str(), (channel_name+" workspace").c_str()); RooWorkspace &proto = *protoOwner; - auto proto_config = make_unique("ModelConfig", &proto); + auto proto_config = std::make_unique("ModelConfig", &proto); // preprocess functions for(auto const& func : fPreprocessFunctions){ @@ -814,8 +813,8 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo std::vector> allSampleHistFuncs; std::vector sampleScaleFactors; - std::vector< pair > statNamePairs; - std::vector< pair> > statHistPairs; // + std::vector> statNamePairs; + std::vector>> statHistPairs; // const std::string statFuncName = "mc_stat_" + channel_name; string prefix; @@ -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 @@ -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 @@ -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 @@ -1210,11 +1206,11 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo // Create the histogram of (binwise) // stat uncertainties: - unique_ptr fracStatError( MakeScaledUncertaintyHist( channel_name + "_StatUncert" + "_RelErr", statHistPairs) ); + std::unique_ptr 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, @@ -1281,9 +1277,9 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo for(unsigned int i=0; i( - ("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( + ("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()); @@ -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 asimov_dataset(AsymptoticCalculator::GenerateAsimovData(*model, observables)); + std::unique_ptr asimov_dataset(AsymptoticCalculator::GenerateAsimovData(*model, observables)); proto.import(*asimov_dataset, RooFit::Rename("asimovData")); } @@ -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 @@ -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 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 pdfMap; - vector models; + std::map pdfMap; + vector models; - RooArgList obsList; - for(unsigned int i = 0; i< ch_names.size(); ++i){ - obsList.add(*static_cast(chs[i]->obj("ModelConfig"))->GetObservables()); - } + RooArgList obsList; + for (unsigned int i = 0; i < ch_names.size(); ++i) { + obsList.add(*static_cast(chs[i]->obj("ModelConfig"))->GetObservables()); + } cxcoutI(HistFactory) <<"full list of observables:\n" << obsList << std::endl; RooArgSet globalObs; @@ -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. @@ -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)); @@ -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 ); @@ -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; @@ -1812,7 +1810,4 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo return std::unique_ptr(ErrorHist); } - -} // namespace RooStats -} // namespace HistFactory - +} // namespace RooStats::HistFactory