P0 value computing with StandardHypostestInv

[diboye]

Dear all,

I created a workspace in Roofit from 2D histogram as in this example here:

#include "RooStats/HistFactory/MakeModelAndMeasurementsFast.h"
#include "TH2D.h"



void ShapeFit_diallo()
{
  //Histos were already made by ShapeFit_makeHists.C. Let's open them up and use them
  TFile* file_signal = new TFile("hists/BRscaled_gaussiansignal_allchannel_mZd30_mH110GeV.root");
  TFile* file_data_bkg = new TFile("hists/Data_bkg_2D_all_channel_3-2GeVBin_interpNegBin.root");
  //TFile* file_data = new TFile("hists/data_allchannel.root");
  //signal
  TH2D* h_sig = (TH2D*)file_signal->Get("Nominal/sig");
  //TH2D* h_sig_scale_down = (TH2D*)file_bkg->Get("Nominal/bkg");
  //TH2D* h_sig_scale_up = (TH2D*)file->Get("Signal_scale_up");
  //background
  TH2D* h_bkg = (TH2D*)file_data_bkg->Get("bkg_all");
  //TH2D* h_bkg_slope_down = (TH2D*)file->Get("Background_slope_down");
  //TH2D* h_bkg_slope_up = (TH2D*)file->Get("Background_slope_up");
  //data
  TH2D* h_obs = (TH2D*)file_data_bkg->Get("data_all");

//--------------------

//create the histfactory model
  RooStats::HistFactory::Measurement meas("ShapeFit", "ShapeFit"); 
  meas.SetOutputFilePrefix("workspaces/ShapeFit/2D");
  meas.SetExportOnly(1);
  meas.SetPOI("mu");

  // this scales the histogram content, which already includes lumi, so set to 1
  meas.SetLumi(1.0);
  meas.SetLumiRelErr(0.02);


//create the parameters for the model

  //first the signal normalization
  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

//create the SR
  RooStats::HistFactory::Channel SR("SR_shape");
  SR.SetData(h_obs);

//add the signal and background samples
  RooStats::HistFactory::Sample sample_S("S_shape");
  sample_S.SetHisto(h_sig);
  sample_S.AddNormFactor(normS);
//add the shape uncertainty for the signal
  //RooStats::HistFactory::HistoSys shapeSys_S("Scale");
  //shapeSys_S.SetHistoHigh(h_sig_scale_up);
  //shapeSys_S.SetHistoLow(h_sig_scale_down);
  //sample_S.AddHistoSys(shapeSys_S);
//let's also add a normalization systematic
  //sample_S.AddOverallSys("sys_S",1./1.05,1.05);

  SR.AddSample(sample_S);

  RooStats::HistFactory::Sample sample_B("B_shape");
  sample_B.SetHisto(h_bkg);
//add the shape uncertainty for the signal
  //RooStats::HistFactory::HistoSys shapeSys_B("Slope");
  //shapeSys_B.SetHistoHigh(h_bkg_slope_up);
  //shapeSys_B.SetHistoLow(h_bkg_slope_down);
  //sample_B.AddHistoSys(shapeSys_B);
//also add a normalization uncertainty
  //sample_B.AddOverallSys("sys_B",1./1.15,1.15);
  //SR.AddSample(sample_B);

//add the single region to the measurement
  meas.AddChannel(SR);

//make the workspace
  RooStats::HistFactory::MakeModelAndMeasurementFast(meas);
}

Then I wanted to compute the p0 value using StandardHypostestInv from this: ROOT: tutorials/roostats/StandardHypoTestInvDemo.C File Reference . I copy the macro and gives it my workspace as input but when I try to run it nothing happen. I attached my workspace.
Could someone tell me please what I am missing?

Thanks

Best,

Diallo.
2D_combined_ShapeFit_model.root (54.3 KB)

[moneta]

Hi,
What you would like to do, compute a limit or performing an hypothesis test ?

Looking at your workspace, I see a problem with your input pdf describing the shape, it seems you have only the signal shape and this histogram is full of zero bins. You are missing the background histogram. Note that the overall pdf cannot be zero in some bins, otherwise you will get an error, when computing the likelihood.

Lorenzo

[diboye]

Hi,

Thanks for replying. What I want to do is to first create a workspace from the 2D histograms (of bkg, signal and data) using the code above. I am also attaching the root file containing the 2D histograms (which are not empty) that I am using as input for the code. Then compute the p0,so the hypothesis test. But now I am even wondering if my workspace is well constructed.
BRscaled_gaussiansignal_allchannel_mZd30_mH110GeV.root (108.2 KB)
Data_bkg_2D_all_channel_3-2GeVBin_interpNegBin.root (10.8 KB)

[moneta]

Hi,

What I was saying that you have included only the signal bin in the workspace, and probably for this is full of bins with zero bin content. You need to include also the background histogram as input to the HistFactory.

Lorenzo

[diboye]

Hi,

That’s right, I was missing that, fixed now. Thanks

[diboye]

Hi Lorenzo,

Does the workspace I created seems right?
I run the StandarHypoTestInv with 10k toys but the result doesn’t seem to be convincing. Here is what I got:

Results HypoTestCalculator_result: 
 - Null p-value = 0 +/- 0
 - Significance = inf +/- -nan sigma
 - Number of Alt toys: 0
 - Number of Null toys: 10000
 - Test statistic evaluated on data: 1.78464
 - CL_b: 0 +/- 0
 - CL_s+b: 0 +/- inf
Error: Cannot compute CLs because CLb = 0. Returning CLs = -1
 - CL_s: -1 +/- -1
Result for workspace file 2D_combined_ShapeFit_model.root mu fitted = 1.12858e-06 +/- 9.60453e-07 null=0 alt=1 in range 0:100 and test stat #2
10000 null samples: 0 0 0 0 0 0 0 0 0 0...
10000 null toys, 0 fit failures
Null p-value 0 +/- 0
Asymptotic p-value 0.0294288, significance 1.88925 (1DoF, 1-sided)

I attached again the workspace after after the bkg.
I am running now for 25k toys to see it will make any difference.

Best,

Diallo.
2D_combined_ShapeFit_model.root (71.6 KB)

[moneta]

Hi,
Your model is still zero in some allowed region of the observables. This will cause all fit to fail.
Here is a plot of your model pdf in the two observables:

[diboye]

Hi,

I guess you are talking about the blank region? Actually the observables are defined in such a way that obs_y (the child particle) <= obs_x/2 (the parent particle) that’s what makes the triangular shape. And Also 10<obs_x<50 and obs_y < 115. So should I add this constraints in macro that created the workspace I guess?