Excluding signal region from fit and background subtraction

[SimN]

Hi,

I’m trying to fit a histogram in two regions while excluding the signal region. and then I would like to subtract the fit function from the original histogram.
The fit is working perfectly but the problem is the background subtraction was not successful.
You can find my code and the pictures attached.
When I fit the full histogram, the background subtraction works. Could you please tell me what I’m doing wrong here?

 RooRealVar x("x","MM_{(#pi #pi)} GeV/c^{2}",3.04,3.17) ;

  RooDataHist dh("dh","dataset",x,h1);
  RooPlot* frame = x.frame(Title("2regions")) ;
  x.setRange("Range1",3.04,3.089);
  x.setRange("Range2",3.11,3.17);
  dh.plotOn(frame,Name("dh")) ;

  RooRealVar a0("a0","a0",-2.,2.) ;
  RooRealVar a1("a1","a1",-2.,2.) ;
  RooRealVar a2("a2","a2",-2.,2.) ;

  RooChebychev bkg("bkg","background p.d.f.",x,RooArgList(a0));
  
  RooRealVar nbkg("N_{BKG}","signal background even0ts",0,100000000);
RooAddPdf all("all","model",RooArgList(bkg),RooArgList(nbkg));
  RooChi2Var chi2("chi2","chi2",all,dh,Range("Range1,Range2"),NumCPU(2));
  RooMinuit m(chi2);
  m.migrad();
  m.hesse();
  RooFitResult* r = m.save();
all.plotOn(frame,Range("Range1,Range2"),NormRange("Range1,Range2"),LineColor(kRed),Components(bkg));
  all.paramOn(frame,Layout(0.5,0.90,0.55));
  TH1 *hbkg = bkg.createHistogram("hbkg",x,Binning(150,3.04,3.17));
  frame->Draw();
c2->Divide(1,2);
c2->cd(1);
hbkg->Scale(nbkg.getVal());
hbkg->Sumw2();
hbkg->Draw();
hbkg->SetTitle("Background histogram");
c2->cd(2);
TH1F *h4260 = (TH1F*)hist4260->Clone("h4260");

h4260->Sumw2();
h4260->Add(hbkg,-1);
h4260->Draw();
h4260->SetTitle("Background subtracted distribution");
h4260->GetXaxis()->SetTitle("MM_{(#pi#pi)} GeV/c^{2}");

TFile *file = new TFile("4260_2RG.root", "RECREATE");
h4260->Write();

                                                             

[amadio]

The normalization of your histograms seems to be off. If you look at the second histogram, it has many more counts than the original. That might be due to different binning. However, in this particular case, if you just want to subtract the background, I recommend using a simple line instead of another histogram for the subtraction.

[SimN]

hi @amadio,
Thank you for your reply. I realized that there’s something wrong about the normalization but I don’t really know exactly how to fix it.
Concerning the background subtraction, I thought this was the only way to do it in Roofit. Fit with a function and then create a histogram using (createHistogram) out of the function.
Could you please give me more details on how to do a simple line subtraction ?

[StephanH]

Hi @SimN,

as @amadio pointed out, the normalisation seems to be off. That could be because nbkg does not refer to the full range of your data, but maybe just to Range2, the one you defined last. For these cases, there is
fixCoefRange()

You can use it to make clear that nbkg refers to the full range of your observable, although you only fit the sidebands. That might yield the number you actually need to correctly do the subtraction.

[SimN]

Hi @StephanH and @amadio,

Don’t know if you could help me with another problem as well.
So after fitting the histograms in two regions and subtracting the background. I get the subtracted histogram but all my attempts to fit it are failing. I tried chi-square fit and likelihood fit, nothing worked.
While if I subtract the full range, I can easily fit the background subtracted histogram using a chi-square fit.
You can see both histograms (and the root files) attached. They look the same, but there’s difference in the number of entries.Any ideas on what I’m doing wrong here?

4230_2RG.root (4.9 KB)
4230_BS.root (4.9 KB)

My code :

 RooRealVar x("x","MM_{(#pi #pi)} GeV/c^{2}",3.04,3.17) ;

  RooDataHist dh("dh","dataset",x,h1);
  RooPlot* frame = x.frame(Title("@4.26 GeV :: e^{+}e^{-} -> J/#psi #pi^{+} #pi^{-}")) ;
  x.setRange("Range1",3.04,3.089);
  x.setRange("Range2",3.11,3.17);
  x.setRange("fullrange",3.04,3.17);


  dh.plotOn(frame,Name("dh")) ;


  RooRealVar a0("a0","a0",-2.,2.) ;
  RooRealVar a1("a1","a1",-2.,2.) ;
  RooRealVar a2("a2","a2",-2.,2.) ;

  RooChebychev bkg("bkg","background p.d.f.",x,RooArgList(a0));
  RooChebychev bkg1("bkg1","background p.d.f.",x,RooArgList(a0,a1));
  RooChebychev bkg2("bkg2","background p.d.f.",x,RooArgList(a0,a1,a2));


  RooRealVar nsig("N_{SIG}","signal events",0,1000000);
  RooRealVar nbkg("N_{BKG}","signal background even0ts",0,100000000);

  RooAddPdf all("all","all",RooArgList(bkg),RooArgList(nbkg));
  RooFitResult* r = all.fitTo(dh,Extended(kTRUE),Save()) ;  all.plotOn(frame,Range("Range1,Range2"),NormRange("Range1,Range2"),LineColor(kRed),Components(bkg));
  all.paramOn(frame,Layout(0.5,0.90,0.55));
 
 TH1 *hbkg = bkg.createHistogram("hbkg",x,Binning(150,3.04,3.17));
  frame->Draw();

c2->Divide(1,2);
c2->cd(1);
hbkg->Scale(nbkg.getVal());
hbkg->Sumw2();
hbkg->Draw();
hbkg->SetTitle("Background histogram");

c2->cd(2);
TH1F *h4230 = (TH1F*)hist4230->Clone("h4230");
h4230->Sumw2();
h4230->Add(hbkg,-1);
h4230->Draw();
h4230->GetXaxis()->SetTitle("MM_{(#pi#pi)} GeV/c^{2}");
TFile *file = new TFile("4230_2RG.root", "RECREATE");
h4230->Write();
file->Close();

[StephanH]

• I don’t know what you mean with subtract full range. There doesn’t seem to be a difference between the histograms except for the number of entries, so it’s not clear which ranges are meant here.

• As to the number of entries, it’s obviously difficult for ROOT to correctly compute the number of entries. We are dealing with the difference of a histogram where entries are simply counted and a histogram where entries are made up from a shape and a normalisation. What is the statistical uncertainty of the entries in each “made up” bin? Do you expect ROOT to correctly compute the entries if you subtract the two histograms? Do you even expect the uncertainties in each bin to be correct?

• Hint: Look for SumW2 in the TH1 documentation. You will definitely need that if you want to proceed on the subtraction path.

• I wanted to say this before, but at the time I thought that maybe it’s more complicated than what you asked for. Now, I think that what I wanted to propose is much easier:
  Don’t subtract, make a combined model. Use a RooAddPdf to sum a signal and a background model, and do a combined fit. This will treat all uncertainties correctly, and you will get a reasonable result.