Skip to content

Commit 1446a30

Browse files
committed
[RF] Don't ignore ParamSetting values in HistFactory Asimov dataset
HistFactory allows to set nuisance parameters to given values (and possibly set them to constant) by using something like ``` <ParamSetting Val=1 Const="True"> myNP </ParamSetting> ``` However the Asimov dataset created by `hist2workspace` does not respect this setting, because in [MakeCombiendModel](https://root.cern/doc/master/HistoToWorkspaceFactoryFast_8cxx_source.html#l02109) the Asimov workspace is created before the ParamSettings are applied. This commit is fixing this problem. The problem can be reproduced with the following script/instructions: ```C++ using namespace RooStats; using namespace HistFactory; // Step 1: initialize the XML files // // Run prepareHistFactory before in this directory // // Modify the GaussExample in the generated example to be: // // <Measurement Name="GaussExample" Lumi="1." LumiRelErr="0.1" > // <POI>SigXsecOverSM</POI> // <ParamSetting Const="True">Lumi</ParamSetting> // <ParamSetting Val="20" Const="True">alpha_syst1</ParamSetting> // </Measurement> // Step 2: Create the RooWorkspace // ConfigParser parser{}; auto measurements = parser.GetMeasurementsFromXML("config/example.xml"); // Get the first example, the GaussExample Measurement& meas = measurements.front(); // Collect the histograms from their files, meas.CollectHistograms(); // Now, do the measurement std::unique_ptr<RooWorkspace> ws{MakeModelAndMeasurementFast( meas )}; // Step 3: Validation plot to show mismatch between model and Asimov data // auto x = ws->var("obs_x_channel1"); auto data = ws->data("asimovData"); auto pdf = ws->pdf("model_channel1"); // Draw the asimov data together with the model and check if they match TCanvas c1("c1"); RooPlot *xframe = x->frame(); data->plotOn(xframe); pdf->plotOn(xframe); xframe->Draw(); c1.SaveAs("plot.png"); // Note that the errors will be plotted wrongly if the Asimov data is // correct. That's because if the changed alpha_syst1 value is considered // correctly, the asimov dataset has non-integer weights. In this case, the // plotting uses SumW2 errors instead of Poisson erros, and the weights // squared is not filled correctly in the Asimov data that is created by // AsymptoticCalculator::GenerateAsimovData(). But this is a completely // different issue. ```
1 parent 565e9aa commit 1446a30

File tree

2 files changed

+30
-33
lines changed

2 files changed

+30
-33
lines changed

roofit/histfactory/inc/RooStats/HistFactory/HistoToWorkspaceFactoryFast.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -100,7 +100,7 @@ namespace RooStats{
100100
std::vector<std::unique_ptr<RooWorkspace>>& wspace_vec,
101101
std::vector<std::string> const& channel_names,
102102
std::string const& dataSetName,
103-
RooArgList obsList,
103+
RooArgList const& obsList,
104104
RooCategory* channelCat);
105105

106106
RooHistFunc* MakeExpectedHistFunc(const TH1* hist, RooWorkspace* proto, std::string prefix,

roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx

Lines changed: 29 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -1872,55 +1872,29 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
18721872
combined->defineSet("globalObservables",globalObs);
18731873
combined_config->SetGlobalObservables(*combined->set("globalObservables"));
18741874

1875-
1876-
////////////////////////////////////////////
1877-
// Make toy simultaneous dataset
1878-
cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1879-
<< "\tcreate toy data for " << channelString.str()
1880-
<< "\n-----------------------------------------\n" << endl;
1881-
1882-
1883-
// now with weighted datasets
1884-
// First Asimov
1885-
//RooDataSet * simData=nullptr;
18861875
combined->factory("weightVar[0,-1e10,1e10]");
18871876
obsList.add(*combined->var("weightVar"));
1877+
combined->defineSet("observables",{obsList, *channelCat}, /*importMissing=*/true);
1878+
combined_config->SetObservables(*combined->set("observables"));
18881879

1889-
// Create Asimov data for the combined dataset
1890-
std::unique_ptr<RooDataSet> asimov_combined{static_cast<RooDataSet*>(AsymptoticCalculator::GenerateAsimovData(*simPdf,
1891-
obsList))};
1892-
if( asimov_combined ) {
1893-
combined->import( *asimov_combined, Rename("asimovData"));
1894-
}
1895-
else {
1896-
std::cout << "Error: Failed to create combined asimov dataset" << std::endl;
1897-
throw hf_exc();
1898-
}
18991880

19001881
// Now merge the observable datasets across the channels
19011882
for(RooAbsData * data : chs[0]->allData()) {
1902-
// Only include RooDataSets where a data with the same name doesn't only
1903-
// exist in the merged workspace (like the "asimovData").
1904-
if(dynamic_cast<RooDataSet*>(data) && !combined->data(data->GetName())) {
1883+
// We are excluding the Asimov data, because it needs to be regenerated
1884+
// later after the parameter values are set.
1885+
if(std::string("asimovData") != data->GetName()) {
19051886
MergeDataSets(combined, chs, ch_names, data->GetName(), obsList, channelCat);
19061887
}
19071888
}
19081889

19091890

1910-
obsList.add(*channelCat);
1911-
combined->defineSet("observables",obsList);
1912-
combined_config->SetObservables(*combined->set("observables"));
1913-
19141891
if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::INFO))
19151892
combined->Print();
19161893

19171894
cxcoutP(HistFactory) << "\n-----------------------------------------\n"
19181895
<< "\tImporting combined model"
19191896
<< "\n-----------------------------------------\n" << endl;
19201897
combined->import(*simPdf,RecycleConflictNodes());
1921-
//combined->import(*simPdf, RenameVariable("SigXsecOverSM","SigXsecOverSM_comb"));
1922-
// cout << "check pointer " << simPdf << endl;
1923-
// cout << "check val " << simPdf->getVal() << endl;
19241898

19251899
std::map< std::string, double>::iterator param_itr = fParamValues.begin();
19261900
for( ; param_itr != fParamValues.end(); ++param_itr ){
@@ -1957,6 +1931,29 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
19571931
combined->importClassCode();
19581932
// combined->writeToFile("results/model_combined.root");
19591933

1934+
1935+
////////////////////////////////////////////
1936+
// Make toy simultaneous dataset
1937+
cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1938+
<< "\tcreate toy data for " << channelString.str()
1939+
<< "\n-----------------------------------------\n" << endl;
1940+
1941+
1942+
// now with weighted datasets
1943+
// First Asimov
1944+
1945+
// Create Asimov data for the combined dataset
1946+
std::unique_ptr<RooDataSet> asimov_combined{static_cast<RooDataSet*>(AsymptoticCalculator::GenerateAsimovData(
1947+
*combined->pdf("simPdf"),
1948+
obsList))};
1949+
if( asimov_combined ) {
1950+
combined->import( *asimov_combined, Rename("asimovData"));
1951+
}
1952+
else {
1953+
std::cout << "Error: Failed to create combined asimov dataset" << std::endl;
1954+
throw hf_exc();
1955+
}
1956+
19601957
return combined;
19611958
}
19621959

@@ -1965,7 +1962,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
19651962
std::vector<std::unique_ptr<RooWorkspace>>& wspace_vec,
19661963
std::vector<std::string> const& channel_names,
19671964
std::string const& dataSetName,
1968-
RooArgList obsList,
1965+
RooArgList const& obsList,
19691966
RooCategory* channelCat) {
19701967

19711968
// Create the total dataset

0 commit comments

Comments
 (0)