Basic questions about Histfactory

[june0812]

Dear expert,

I recently started to work with Histfactory. I try to write simple code and test the functions of Histfactory. My code to make a workspace is

	// data
	TH1F* data_hist = new TH1F("data_hist", "data_hist", 2, 0, 2);
	for (int i = 0; i < 10000; i++) data_hist->Fill(0.5,0.01);
	for (int i = 0; i < 1000; i++) data_hist->Fill(1.5, 0.1);

	// MC signal pdf
	TH1F* Signal = new TH1F("Signal", "Signal", 2, 0, 2);
	Signal->SetBinContent(1, 10);
	Signal->SetBinContent(2, 80);

	TH1F* Signal_up = new TH1F("Signal_up", "Signal_up", 2, 0, 2);
	Signal_up->SetBinContent(1, 12);
	Signal_up->SetBinContent(2, 84);
	Signal_up->Add(Signal, -1.0);

	TH1F* Signal_down = new TH1F("Signal_down", "Signal_down", 2, 0, 2);
	Signal_down->SetBinContent(1, 7);
	Signal_down->SetBinContent(2, 81);
	Signal_down->Add(Signal, -1.0);

	// MC Background pdf
	TH1F* Background1 = new TH1F("Background1", "Background1", 2, 0, 2);
	Background1->SetBinContent(1, 20);
	Background1->SetBinContent(2, 20);


	// work space

	// set data
	HistFactory::Data data;
	data.SetHisto(data_hist);

	// add measurement
	HistFactory::Measurement meas("MyAnalysis");
	meas.SetPOI("mu");
    meas.SetLumi(1.0);
    meas.AddConstantParam("Lumi");

	// add channel
	HistFactory::Channel channel("channel");
	channel.SetData(data);

	// add signal sample in channel
	HistFactory::Sample signal_sample("signal_sample");
	signal_sample.SetHisto(Signal);

	RooStats::HistFactory::HistoSys shapeSys_S("shape_signal");
	shapeSys_S.SetHistoHigh(Signal_up);
	shapeSys_S.SetHistoLow(Signal_down);
	signal_sample.AddHistoSys(shapeSys_S);

	RooStats::HistFactory::NormFactor normS;
	normS.SetName("mu");
	normS.SetHigh(100); // maximum value it can take
	normS.SetLow(0); // minimum value it can take
	normS.SetVal(1); // nominal value
	signal_sample.AddNormFactor(normS);

    //  signal_sample.ActivateStatError();
	channel.AddSample(signal_sample);

	// add background sample in chennel
	HistFactory::Sample background_sample("background_sample");
	background_sample.SetHisto(Background1);

    //  background_sample.ActivateStatError();
	channel.AddSample(background_sample);

	// add channel in measurement
	meas.AddChannel(channel);

	HistFactory::HistoToWorkspaceFactoryFast h2w(meas);
	RooWorkspace* w = h2w.MakeCombinedModel(meas);
	w->Print("t");
	w->SetName("w");
	w->writeToFile("test_file.root");

data_hist is my (weighted) data. I want to fit, using Signal and Background1. I also add shapeSys_S to add systematic uncertainty on shape of Signal. My code to fit is

	TFile* f = TFile::Open("test_file.root");
	RooWorkspace* w = (RooWorkspace*)f->Get("w");
	w->Print();
	RooAbsData* data = w->data("obsData");

	ModelConfig* sbModel = (ModelConfig*)w->obj("ModelConfig");
	RooFitResult* fitres = sbModel->GetPdf()->fitTo(*data);

	RooSimultaneous* simPdf = (RooSimultaneous*)sbModel->GetPdf();

    TCanvas* cdata = new TCanvas("sPlot","sPlot demo", 700, 700);
    RooPlot* frame = w->var("obs_x_channel")->frame();
    data->plotOn(frame, RooFit::DataError(RooAbsData::SumW2));
    w->pdf("channel_model")->plotOn(frame);
    frame->Draw();

This code shows

I have several questions

1. The above plot shows strange uncertainty. Because I use

	TH1F* data_hist = new TH1F("data_hist", "data_hist", 2, 0, 2);
	for (int i = 0; i < 10000; i++) data_hist->Fill(0.5,0.01);
	for (int i = 0; i < 1000; i++) data_hist->Fill(1.5, 0.1);

statistical uncertainty should be smaller than the above plot’s one. What should I do?

2. I want to show signal and background components separately. Example plot is
   

   
   
   (source: Measurement of inclusive and differential cross sections in the $H \rightarrow ZZ^* \rightarrow 4\ell$ decay channel in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector)
   How can I draw this kind of plot?

3. If I want to add statistical fluctuation of templates, is it OK to add

signal_sample.ActivateStatError();

and

background_sample.ActivateStatError();

in my first code?

Sorry for basic question. I cannot find detail example code in histfactory manual (https://cds.cern.ch/record/1456844/files/CERN-OPEN-2012-016.pdf). It only talks about how histfactory works. Also, there is no enough tutorial code in root (ROOT: HistFactory Tutorials).

Thank you

[dastudillo]

For 2: if all are histograms, you can add them to a THStack and draw with the option “nostack”; also add them to a TLegend. Instead of a THStack, you can draw them one by one on the same canvas with (for histograms) or without (for TGraphs) the option “same”, as shown in the TLegend doc above.

[june0812]

Sorry for the ambiguous question.

My question is how to draw histogram pdf, not histogram. For histogram pdf, there is not Draw member function.

Alternatively, I found components argument for plotOn function. Here is my fixed code:

  TFile* f = TFile::Open("test_file.root");
  RooWorkspace* w = (RooWorkspace*)f->Get("w");
  w->Print();
  RooAbsData* data = w->data("obsData");

  ModelConfig* sbModel = (ModelConfig*)w->obj("ModelConfig");
  RooFitResult* fitres = sbModel->GetPdf()->fitTo(*data);

  RooSimultaneous* simPdf = (RooSimultaneous*)sbModel->GetPdf();

  TCanvas* cdata = new TCanvas("sPlot", "sPlot demo", 700, 700);
  RooPlot* frame = w->var("obs_x_channel")->frame();
  data->plotOn(frame, RooFit::DataError(RooAbsData::SumW2));
  w->pdf("channel_model")->plotOn(frame, DrawOption("F"), FillColor(kRed));
  w->pdf("channel_model")->plotOn(frame, DrawOption("F"), FillColor(kBlue), Components("L_x_signal_sample_channel_overallSyst_x_HistSyst"));
  frame->Draw();

However, I found some problem.
For ROOT Version: 6.16/00, the above fixed code works fine:

red histogram is Background1, and blue histogram is Signal. Though, when I run the same code at ROOT Version: 6.24/00, It shows

blue and red histograms are normalized by the number of data, respectively. I carefully guess it is a bug, because there was a similar bug at past (HistFactory plotting of components broken in ROOT 6.16/00).