Dear Rooters,
I have a problem in the plotting of the results of a fit that to me seems a bug introduces in only when using the ROOT version 6.14.
We are fitting a 3D data histogram using a model constructed with the RooHistFactory package, and at the end of the fit we are plotting the results using the RooAbsPdf::plotOn() method. In order to get a histogram with all the fit components stacked one on top of each other, we are looping over a set of component names and using the command RooFit::Components() to draw only a subset of the components.
We do this recursevely to plot all of them. To be more clear, if we have three components we do
model->plotOn(frame, Components("first"), MoveToBack());
model->plotOn(frame, Components("first,second"), MoveToBack());
model->plotOn(frame, Components("first,second,third"), MoveToBack());
As long as we use ROOT v6.08, the plot looks fine and all the components are picked up correctly.
If we run the same code, though, with the ROOT 6.14, the code does not complain about components not found, but each component is drawn (both in normalization and shape) as the full model.
To reproduce the error, I have modified one of the histfactory tutorials in order to include most of the features we are including in our analysis (in case it can help to debug).
#include "TTree.h"
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "RooPolynomial.h"
#include "RooCustomizer.h"
#include "RooChi2Var.h"
#include "RooAbsData.h"
#include "RooRealSumPdf.h"
#include "RooPoisson.h"
#include "RooGaussian.h"
#include "RooRealVar.h"
#include "RooMCStudy.h"
#include "RooMinuit.h"
#include "RooCategory.h"
#include "RooHistPdf.h"
#include "RooSimultaneous.h"
#include "RooExtendPdf.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooFitResult.h"
#include "RooMsgService.h"
#include "RooParamHistFunc.h"
#include "RooHist.h"
#include "RooRandom.h"
#include "RooPlot.h"
#include "RooCmdArg.h"
#include "TString.h"
#include "TH1D.h"
#include "RooStats/HistFactory/FlexibleInterpVar.h"
#include "RooStats/HistFactory/PiecewiseInterpolation.h"
#include "RooStats/HistFactory/HistFactorySimultaneous.h"
#include "RooStats/HistFactory/Channel.h"
#include "RooStats/HistFactory/MakeModelAndMeasurementsFast.h"
#include "RooStats/HistFactory/Measurement.h"
#include "RooStats/HistFactory/ParamHistFunc.h"
#include "RooStats/HistFactory/HistFactoryModelUtils.h"
#include "RooStats/HistFactory/RooBarlowBeestonLL.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/MinNLLTestStat.h"
#include <boost/program_options/option.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/variables_map.hpp>
#include <fstream>
#include "RooTFnBinding.h"
#include "RooUniform.h"
#include "TSystem.h"
using namespace RooStats;
using namespace HistFactory;
/*
A ROOT script demonstrating
an example of writing a HistFactory
model using c++ only.
This example was written to match
the example.xml analysis in
$ROOTSYS/tutorials/histfactory/
Written by George Lewis, modified (Simone Meloni)
*/
//void plotProjection(RooCategory* idx,RooAbsData* data,RooAbsPdf* model,RooRealVar* obs,TString folder, TString fitvar,RooStats::ModelConfig* mc,RooWorkspace* w, TString output_file_name, TFile* fin, const char* name_suffix)
void plotProjection(RooCategory* idx, RooAbsData* data, RooSimultaneous* model, RooRealVar* obs)
{
RooPlot *frame = new RooPlot();
frame->SetName(TString(obs->GetName())+TString("_frame_")+idx->getLabel());
frame = obs->frame("");
RooPlot *pframe = new RooPlot();
pframe->SetName(TString(obs->GetName())+TString("_pframe_")+idx->getLabel());
pframe = obs->frame(""); //frame used to plot the pulls
//Define the colour of the component
std::vector<int> colours;
colours.push_back(4);
colours.push_back(38);
colours.push_back(kMagenta);
//Define the names of the components
std::vector<const char*> component_names;
component_names.push_back("hNDpMuNu_");
component_names.push_back("hNDstTauNu_");
component_names.push_back("hNDpTauNu_");
//Try to see all the components (why shoudlnt you?)
TString components = "";
data->plotOn(frame,RooFit::MarkerSize(0.4), RooFit::DrawOption("ZP"), RooFit::DataError( RooAbsData::Poisson), RooFit::Cut(TString("channelCat==channelCat::")+TString(idx->getLabel())));
int colour_count=0;
for (int i = 0; i < component_names.size(); ++i)
{
std::cout << component_names[i] << std::endl;
if (components != "")
{
components+=",";
}
components+="*";
components+=component_names[i];
components+="*";
std::cout << components << " " << colours[colour_count] << std::endl;
model->plotOn(frame, RooFit::Slice(*idx), RooFit::ProjWData(*idx,*data), RooFit::DrawOption("F"),RooFit::LineColor(colours[colour_count]),RooFit::FillColor(colours[colour_count]),RooFit::Components(components),RooFit::MoveToBack(), RooFit::Name(component_names[i]));
colour_count++;
}
//-----Construct the pulls plot
model->plotOn(frame, RooFit::Slice(*idx), RooFit::ProjWData(*idx,*data),RooFit::DrawOption("L"),RooFit::LineColor(kBlue), RooFit::LineWidth(0.000001));
RooHist *pulls = frame->pullHist();
pulls->SetFillColor(kGray+1);
pulls->SetLineColor(kWhite);
pulls->SetMarkerSize(0.01);
pframe->addPlotable(pulls, "B");
//-----------Draw the fit result with the pulls plot below
TCanvas *c1 = new TCanvas("c1", "c1");
c1->Divide(1,2);
//----Dimensions of the single pads
double xmin_1 = 0.005;
double xmax_1 = 0.995;
double ymin_1 = 0.35;
double ymax_1 = 0.995;
double xmin_2 = .005;
double xmax_2 = .995;
double ymin_2 = .05;
double ymax_2 = .35;
c1->cd(1)->SetPad(xmin_1, ymin_1, xmax_1, ymax_1 );
c1->cd(1)->SetBottomMargin(0);
c1->cd(1)->SetTopMargin(0.1);
c1->cd(1)->SetLeftMargin(0.14);
c1->cd(2)->SetPad(xmin_2, ymin_2, xmax_2, ymax_2 );
c1->cd(2)->SetTopMargin(0);
c1->cd(2)->SetBottomMargin(0.14/(ymax_2-ymin_2));
c1->cd(2)->SetLeftMargin(0.14);
//-draw the bottom pad
c1->cd(2);
TPad *padbottom = new TPad("pullspad", "pullspad", 0., 0., 1., 0.3);
pframe->SetTitle("");
pframe->GetYaxis()->SetTitle("Pulls");
pframe->GetYaxis()->SetNdivisions(5);
pframe->GetXaxis()->SetNdivisions(frame->GetXaxis()->GetNdivisions());
pframe->GetYaxis()->SetTitleOffset(0.95*(ymax_2-ymin_2));
pframe->GetYaxis()->SetTitleSize(0.060/(ymax_2-ymin_2));//rescale the title and the label sizes after dividing the pad
pframe->GetYaxis()->SetLabelSize(0.060/(ymax_2-ymin_2));
pframe->GetXaxis()->SetLabelSize(0.060/(ymax_2-ymin_2));
pframe->GetXaxis()->SetTitleSize(0.060/(ymax_2-ymin_2));
pframe->SetMaximum(7.);
pframe->SetMinimum(-7.);
Double_t xmin = pulls->GetMinimum();
Double_t xmax = pulls->GetMaximum();
TF1* upperLine_function = new TF1("upperLine", "3.", xmin, xmax);
TF1* mediumLine_function= new TF1("mediumLine", "0.", xmin, xmax);
TF1* lowerLine_function = new TF1("lowerLine", "-3.",xmin, xmax);
RooAbsReal* upperLine_roofit = RooFit::bindFunction(upperLine_function, *obs);
RooAbsReal* mediumLine_roofit= RooFit::bindFunction(mediumLine_function,*obs);
RooAbsReal* lowerLine_roofit = RooFit::bindFunction(lowerLine_function, *obs);
upperLine_roofit->plotOn(pframe, RooFit::LineStyle(3), RooFit::LineColor(12));
mediumLine_roofit->plotOn(pframe, RooFit::LineStyle(3),RooFit::LineColor(12));
lowerLine_roofit->plotOn(pframe, RooFit::LineStyle(3), RooFit::LineColor(12));
pframe->Draw();
//--draw the top pad
c1->cd(1);
frame->GetYaxis()->SetTitleSize(0.060/(ymax_1-ymin_1)); //rescale the title and the label sizes after dividing the pad
frame->GetYaxis()->SetLabelSize(0.060/(ymax_1-ymin_1));
frame->GetYaxis()->SetTitleOffset(0.95*(ymax_1-ymin_1));
frame->SetMinimum(0.01);
//c1->cd(1)->SetLogy();
frame->Draw();
c1->Print(TString("fit_result")+TString("_")+TString(idx->getLabel())+TString(".pdf"));
}
int main() {
//std::string InputFile = "./data/example.root";
std::string InputFile = "./data/ShapeSys2D.root";
// in case the file is not found
bool bfile = gSystem->AccessPathName(InputFile.c_str());
if (bfile) {
std::cout << "Input file is not found - run prepareHistFactory script " << std::endl;
gROOT->ProcessLine(".! prepareHistFactory .");
bfile = gSystem->AccessPathName(InputFile.c_str());
if (bfile) {
std::cout << "Still no " << InputFile << ", giving up.\n";
exit(1);
}
}
// Create the measurement
Measurement meas("meas", "meas");
meas.SetOutputFilePrefix( "./results/example_UsingC" );
meas.SetPOI( "SigXsecOverSM" );
meas.AddConstantParam("alpha_syst1");
meas.AddConstantParam("Lumi");
meas.SetLumi( 1.0 );
meas.SetLumiRelErr( 0.10 );
meas.SetExportOnly( kTRUE );
//meas.SetBinHigh( 2 );
// Create a channel
Channel chan( "channel1" );
chan.SetData( "data", InputFile );
chan.SetStatErrorConfig( 0.05, "Poisson" );
// Now, create some samples
// Create the signal sample
Sample *signal = new Sample( "hNDpTauNu_sample", "signal", InputFile );
//signal->SetName("hNDpTauNu_channel1");
signal->AddOverallSys( "syst1", 0.95, 1.05 );
signal->AddNormFactor( "SigXsecOverSM", 1, 0, 3 );
chan.AddSample( *signal );
// Background 1
Sample *background1 = new Sample( "hNDpMuNu_sample", "background1", InputFile );
//background1->SetName("hNDpMuNu_channel1");
background1->ActivateStatError();
background1->AddOverallSys( "syst2", 0.95, 1.05 );
chan.AddSample( *background1 );
// Background 1
Sample *background2 = new Sample( "hNDstTauNu_sample", "background2", InputFile );
//background2->SetName("hNDstTauNu_channel1");
background2->ActivateStatError();
background2->AddOverallSys( "syst3", 0.95, 1.05 );
chan.AddSample( *background2 );
// Done with this channel
// Add it to the measurement:
meas.AddChannel( chan );
Channel chan2 (chan);
chan2.SetName("channel2");
meas.AddChannel(chan2);
Channel chan3 (chan);
chan3.SetName("channel3");
meas.AddChannel(chan3);
Channel chan4 (chan);
chan4.SetName("channel4");
meas.AddChannel(chan4);
Channel chan5 (chan);
chan5.SetName("channel5");
meas.AddChannel(chan5);
// Collect the histograms from their files,
// print some output,
meas.CollectHistograms();
meas.PrintTree();
meas.SetExportOnly(kTRUE);
// One can print XML code to an
// output directory:
// meas.PrintXML( "xmlFromCCode", meas.GetOutputFilePrefix() );
// Now, do the measurement
//###### THIS IS THE PART I MODIFIED FOR PLOTTING#######
RooWorkspace *w;
w = MakeModelAndMeasurementFast( meas );
RooStats::ModelConfig *mc = (RooStats::ModelConfig*) w->obj("ModelConfig");
RooArgSet *obs = (RooArgSet*) mc->GetObservables();
RooCategory *idx = (RooCategory*) obs->find("channelCat");
RooRealVar *x1 = (RooRealVar*) obs->find("obs_x_channel1");
RooRealVar *x2 = (RooRealVar*) obs->find("obs_x_channel2");
RooRealVar *x3 = (RooRealVar*) obs->find("obs_x_channel3");
RooRealVar *x4 = (RooRealVar*) obs->find("obs_x_channel4");
RooRealVar *x5 = (RooRealVar*) obs->find("obs_x_channel5");
RooAbsData *data = (RooAbsData*) w->data("obsData");
RooSimultaneous *model = (RooSimultaneous*)mc->GetPdf();
RooSimultaneous *model_ = new RooSimultaneous("simPdf_modified","simPdf_modified", *idx);
RooCustomizer* cust_;
RooAbsPdf* pdf_;
RooRealVar *syst2 = (RooRealVar*) mc->GetNuisanceParameters()->find("alpha_syst2");
RooRealVar *newsyst2 = new RooRealVar("alpha_newsyst2", "alpha_newsyst2", 1., 0.80,1.20);
for (unsigned int j = 0; j < 5; j++) {
idx->setIndex(j);
pdf_ = model->getPdf(idx->getLabel());
cust_ = new RooCustomizer(*pdf_,"cust_");
cust_->replaceArg(*syst2,*newsyst2);
model_->addPdf((RooAbsPdf&)*cust_->build(),idx->getLabel());
}
HistFactorySimultaneous* model_hf = new HistFactorySimultaneous(*model_);
RooAbsReal* nll_hf = model_hf->createNLL(*data,RooFit::Offset(kTRUE));
RooMinuit* minuit_hf = new RooMinuit(*nll_hf);
minuit_hf->setErrorLevel(0.5);
minuit_hf->setStrategy(2);
minuit_hf->fit("smh");
for (unsigned int i = 0; i< 5; i++){
idx->setIndex(i);
if (i==0){
plotProjection(idx, data, (RooSimultaneous*) model_, x1 );
}
if (i==1){
plotProjection(idx, data, (RooSimultaneous*) model_, x2 );
}
if (i==2){
plotProjection(idx, data, (RooSimultaneous*) model_, x3 );
}
if (i==3){
plotProjection(idx, data, (RooSimultaneous*) model_, x4 );
}
if (i==4){
plotProjection(idx, data, (RooSimultaneous*) model_, x5 );
}
}
//model->Print("t");
return 0;
}
which can be copied in the (I assumed under the name example.C) $ROOTSYS/tutorials/histfactory folder and compiled and run with the following python script.
import sys
import os
import ROOT
from ROOT import *
os.system("g++ -o example example.C `root-config --cflags --libs` -lstdc++ -lMinuit -lFoam -lRooFit -lRooFitCore -lRooStats -lHistFactory -lMathMore")
os.system("./example")
The code outputs the following figure when run with ROOT 6.08 (which displays the three fit components correctly)
fit_result_channel1.pdf (14.2 KB)
When run with ROOT 6.14 instead it outputs the following one, in which only the full model with the color assigned to the last component is drawn
fit_result_channel1_6_14.pdf (14.1 KB)
One bit of evidence for a potential bug is that this only happens when working with 2D/3D models in Roofit, the plot looks fine when using a 1D model.
Lastly the code gives a very bad error on exit that starts with
*** glibc detected *** ./example: double free or corruption (out): 0x0000000003b7a760 ***
and ends with a very long backtrace.
Can I ask for any help?
Thanks in advance,
Simone
