Limit setting with HistFactory: convert 2d hist to 1D and use it as input

Dear experts,

I would like to use HistFactory to set limits. In our analysis, we are looking let’s a parent particle Y decaying to two child particles X which decay two 2l each: Y->XX->4l. One way we think we could proceed is to produce first 2D histograms (mll vs m4l, where mll is the average di-lepton mass of the two X particles and m4l is the mass of the parent particle Y). and then convert the 2D histograms to 1D histograms in order to use them as inputs. The issue is when I run the “runAsymptoticsCLs.C” code it never finish. Without throwing any error I am struggling to detect the issue. Does anyone have any advice? I put one mass point model (where mX = 20 GeV and m4l = 175 GeV) in my cern public area: /afs/cern.ch/user/d/diboye/public.

Thanks in advance.

You should clarify what you mean by “converting TH2 to TH1s”.

It looks like an infinite loop. We would need the macro to understand what’s wrong.

Hi,

Thanks for replying. Here the code I implemented which convert a 2D Hist to a 1D hist.

/************************************************
Script from diallo.boye@cern.h
This script converts a 2D histogram to a 1D histogram.
It uses root files inputs containing 2D histogram in TDirectories 
and the outputs contains 1D histogram in the same TDirectories.
Sepember 30th.
It's run by doing: ./command_to_run.sh.
*/
#include "TROOT.h"    
#include "TObject.h" 
#include "TFile.h"
#include "TTree.h"
#include "TH2.h"
#include "TMath.h"
#include "TGraph.h"
#include "TColor.h"
#include "TCanvas.h"
#include "TLegend.h"
#include <TStyle.h>
#include <TString.h>
#include <iostream>
#include <fstream>
#include <sstream>


using namespace std;
void FuncToGet1DFrom2DHist(int binX, int binY,  TH1D* hist1D, TH2D* hist2D, int Ser)
{
  for(int ix =0; ix< binX; ix++)
    {
      for(int iy =0; iy< binY; iy++)
	{
	  //Ser = ix+(iy*ix);
	  Ser = (ix*binY + iy);
	  hist1D->Fill(Ser+1, hist2D->GetBinContent(ix+1,iy+1));
	}
    }
}

void FuncToGet1DFrom2DHist_up_down(int binX, int binY,  TH1D* hist1D_up, TH1D* hist1D_down, TH2D* hist2D_up,TH2D* hist2D_down, int Ser)
{
  for(int ix =0; ix< binX; ix++)
    {
      for(int iy =0; iy< binY; iy++)
	{
	  //Ser = ix+(iy*ix);
	  Ser = (ix*binY + iy);
	  hist1D_up->Fill(Ser+1, hist2D_up->GetBinContent(ix+1,iy+1));
	  hist1D_down->Fill(Ser+1, hist2D_down->GetBinContent(ix+1,iy+1));
	}
    }
}

void DeletingExistingHist()
{
  TObject *obj;
  TKey *key;
  TIter next(gDirectory->GetListOfKeys());
  while ((key = (TKey *) next()))
    {
      obj = gDirectory->Get(key->GetName());
      key->Delete();
    }
}


void convert2Dto1D()
{
  ifstream myfileSig ("listSig.dat");
  ifstream myfileBkg ("listBkg.dat");
  ifstream myfileData ("listData.dat");
  string line;
  vector<string> lineCointainerSig;
  vector<string> lineCointainerBkg;
  vector<string> lineCointainerData;
  TFile* SigFile[50];
  TFile* BkgFile[5];
  TFile* DataFile[5];
  TH2D *tempHist2D;
  TH2D *tempHist2D_up;
  TH2D *tempHist2D_down;
  int Nbinx;
  int Nbiny;
  int Ser;
  TH1D* tempHist1D = new TH1D("tempHist1D","" , 4600, 0, 46000);
  TH1D* tempHist1D_up = new TH1D("tempHist1D_up","" , 4600, 0, 46000);
  TH1D* tempHist1D_down = new TH1D("tempHist1D_down","" , 4600, 0, 46000);
  const char *Syst[] =
    {
      "Nominal",
      "STAT",
      "LUMI",
      "QCDSCALE",
      "PDF",
      "EL_EFF_ID_TOTAL",   
      "EL_EFF_ISO_TOTAL",	
      "EL_EFF_RECO_TOTAL",	
      "MUON_EFF_ISO_STAT",	
      "MUON_EFF_ISO_SYS",  
      "MUON_EFF_RECO_STAT", 
      "MUON_EFF_RECO_STAT_LOWPT",
      "MUON_EFF_RECO_SYS",	
      "MUON_EFF_RECO_SYS_LOWPT",	
      "MUON_EFF_TTVA_STAT",	
      "MUON_EFF_TTVA_SYS",	
      "PRW_DATASF",	
      "EG_RESOLUTION_ALL",	
      "EG_SCALE_ALL",	
      "MUONS_ID",	
      "MUONS_MS",	
      "MUONS_SAGITTA_RESBIAS",     
      "MUONS_SAGITTA_RHO",	
      "MUONS_SCALE"
    };
  
  //storing signal files
    if (myfileSig.is_open())
    {
      while (getline (myfileSig,line))
	{
	  lineCointainerSig.push_back(line);
	}
      myfileSig.close();
    }
  //storing bkg files
  if (myfileBkg.is_open())
    {
      while (getline (myfileBkg,line))
	{
	  lineCointainerBkg.push_back(line);
	}
      myfileBkg.close();
    }

   //storing data files
  if (myfileData.is_open())
    {
      while (getline (myfileData,line))
	{
	  lineCointainerData.push_back(line);
	}
      myfileData.close();
    }
 
  //converting signal
  
    for (int iSig = 0; iSig < 30; iSig++)
    {
      SigFile[iSig] = new TFile(lineCointainerSig[iSig].c_str(), "update");
      for (int j = 0; j < 24; j++)
	{
	  SigFile[iSig]->cd(TString::Format("%s" ,Syst[j]));
	  if (j <2)
	    {
	      tempHist2D = (TH2D*)gDirectory->Get("sig");
	      Nbinx = tempHist2D->GetNbinsX();
	      Nbiny = tempHist2D->GetNbinsY();
	      FuncToGet1DFrom2DHist(Nbinx, Nbiny,  tempHist1D, tempHist2D, Ser);
	      gDirectory->pwd();
	      DeletingExistingHist();
	      tempHist1D->SetName("sig");
	      tempHist1D->Write("");
	      tempHist1D->Reset("");
	    }
	  else
	    {
	      tempHist2D_up = (TH2D*)gDirectory->Get("sig_up");
	      tempHist2D_down = (TH2D*)gDirectory->Get("sig_down");
	      Nbinx = tempHist2D_up->GetNbinsX();
	      Nbiny = tempHist2D_up->GetNbinsY();
	      FuncToGet1DFrom2DHist_up_down(Nbinx, Nbiny, tempHist1D_up, tempHist1D_down, tempHist2D_up, tempHist2D_down, Ser);
	      gDirectory->pwd();
	      DeletingExistingHist();
	      tempHist1D_up->SetName("sig_up");
	      tempHist1D_down->SetName("sig_down");
	      tempHist1D_up->Write("");
	      tempHist1D_down->Write("");
	      tempHist1D_up->Reset("");
	      tempHist1D_down->Reset("");
	    }
	}
       SigFile[iSig]->Close();
    }
  
  // converting bkg
  
   for (int iBkg = 0; iBkg < 3; iBkg++)
    {
      BkgFile[iBkg] = new TFile(lineCointainerBkg[iBkg].c_str(), "update");
      for (int j = 0; j < 24; j++)
	{
	  BkgFile[iBkg]->cd(TString::Format("%s" ,Syst[j]));
	  if (j <2)
	    {
	      tempHist2D = (TH2D*)gDirectory->Get("bkg");
	      Nbinx = tempHist2D->GetNbinsX();
	      Nbiny = tempHist2D->GetNbinsY();
	      FuncToGet1DFrom2DHist(Nbinx, Nbiny,  tempHist1D, tempHist2D, Ser);
	      gDirectory->pwd();
	      DeletingExistingHist();
	      tempHist1D->SetName("bkg");
	      tempHist1D->Write("");
	      tempHist1D->Reset("");
	      tempHist2D->Reset("");
	    }
	  else
	    {
	      tempHist2D_up = (TH2D*)gDirectory->Get("bkg_up");
	      tempHist2D_down = (TH2D*)gDirectory->Get("bkg_down");
	      Nbinx = tempHist2D_up->GetNbinsX();
	      Nbiny = tempHist2D_up->GetNbinsY();
	      FuncToGet1DFrom2DHist_up_down(Nbinx, Nbiny, tempHist1D_up, tempHist1D_down, tempHist2D_up, tempHist2D_down, Ser);
	      gDirectory->pwd();
	      DeletingExistingHist();
	      tempHist1D_up->SetName("bkg_up");
	      tempHist1D_down->SetName("bkg_down");
	      tempHist1D_up->Write("");
	      tempHist1D_down->Write("");
	      tempHist1D_up->Reset("");
	      tempHist1D_down->Reset("");
	      tempHist2D_up->Reset("");
	      tempHist2D_down->Reset("");
	    }
	}
       BkgFile[iBkg]->Close();
    }

    // converting Data
  
   for (int iData = 0; iData < 3; iData++)
     {
       DataFile[iData] = new TFile(lineCointainerData[iData].c_str(), "update");
       tempHist2D = (TH2D*)DataFile[iData]->Get("Data");
       Nbinx = tempHist2D->GetNbinsX();
       Nbiny = tempHist2D->GetNbinsY();
       FuncToGet1DFrom2DHist(Nbinx, Nbiny,  tempHist1D, tempHist2D, Ser);
       gDirectory->pwd();
       DeletingExistingHist();
       tempHist1D->SetName("Data");
       tempHist1D->Write("");
       tempHist1D->Reset("");
       DataFile[iData]->Close();
     }
}



  /*for (int i =0; i < lineCointainerSig.size(); i++)
  {
    cout << " lines : " << lineCointainerSig[i] << '\n';
  }*/
  //cout << " size " << lineCointainerSig.size() << '\n';

it works fine. I also attached a link of a zip file which contains the code and the input in case you want to run it for checking: https://www.dropbox.com/s/myx7g549w2h8tog/Convert2DTo1DHist.zip?dl=0 you can run it by doing: ./command_to_run.sh.
This code creates the input that I am using to set the limit using HistFactory.
I put the package to set the limit in my git repository with a readme file: https://github.com/dialloboye/Statistical-tool.
Just to mention that I have been using the package to set limits where only mll distribution (a 1D histogram) was considered as input with less binning and it was working fine. So I don’t think the package itself has a problem.

Thanks for the extra details. @moneta may have an idea.

1 Like

Hi,

Without throwing any error I am struggling to detect the issue

You can execute the program inside a debugger, e.g. gdb, stop it with ctrl-C and inspect what it’s doing with the backtrace command. That should tell you where your program gets stuck.

Feel free to ping us again or open a bug report at Issues · root-project/root · GitHub if you think the problem is in ROOT!

Cheers,
Enrico

P.S.
I downloaded the code from dropbox and ran bash command_to_run.sh and it looks like it runs to completion without issues. I am unclear how to proceed to reproduce the problem.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.