Subtraction of two roodatahists

Hello all,

I have two RooDataHists. One of them was created using

RooDataHist dh("dh","e-e+",x,Import(*h1)) ;

where h1 is a pointer to a TH1F type histogram.

TH1F *h1 = new TH1F("h1","e-e+ spectra;Mass;Events",150,0.,600.);

The other is created using a user pdf

RooDataHist *back = bg_bern.generateBinned(x,30388) ;

What I wish to do is subtract “back” from “dh”. (Essentially, dh is the entire spectra, and back is the background)

How do I subtract two RooDataHists? For example, for TH1 I could do something like back -> Add(dh,-1).

Hi again,

you cannot do this directly. You can, however, create a standard ROOT histogram from a PDF:
https://root.cern.ch/doc/master/classRooAbsReal.html#a6659d2c301e5cd65b83ee8c9422c2553

The options for the generation of the histogram are also listed in the documentation I linked above.
When you have the TH1x, subtraction of the histograms should be easy.

Hello @StephanH,

In my program, the two histos that I have are h1, h2.

TH1F *h1 = new TH1F("h1","e-e+ spectra;Mass;Events",200,0.,500.);

Then, I created a histo from the roopdf (bg_bern is my fit pdf)


TH1 *h2 = bg_bern.createHistogram("back", x, Binning(7039,0.,500.));

h1 is the spectrum and h2 is the background, so what I want to do is h1-h2 and plot the resultant histo

Is that

h2 -> Add(h1,-1);
h2 -> Draw();

What I should get as output is just a BW ( or a CB or a Gaus) in the middle. However I get back this (4th panel).


which is h2 itself. Where am I going wrong?
I am new to this, could you please also explain, how do I get the bin size same in both histos.?

In principle, you should have the same binning for both histograms, (e.g. 200, 0., 500.), otherwise TH1::Add may fail.
Then, it seems to me that you want to have “spectrum - background”, so:

TH1F *h12 = new TH1F(*h1);
h12->SetNameTitle("h12", "spectrum - background");
h12->Add(h2, -1.0);
h12->Draw();

with

TH1F *h1 = new TH1F("h1","e-e+ spectra;Mass;Events",200,0.,500.); //7039 entries in total
TH1 *h2 = bg_bern.createHistogram("background;mass;count", x, Binning(200,0.,500.));
         // back -> plotOn(frame4);
      	 TH1F *h12 = new TH1F(*h1);
      	 h12 ->SetNameTitle("h12","signal only");
      	 h12 -> Add(h2,-1.0);

I get this

while h2 looks like this (3rd Panel)

Look at the “h1” bins’ contents … the maximum is around 100 … and then compare them with your “h2” … where the maximum is around 0.025 (so, if you subtract them you will see no real difference).

Does the “Binning” in

TH1 *h2 = bg_bern.createHistogram("back", x, Binning(7039,0.,500.));

create 7039 entries spread across 0 to 500?. The documentation says, it’s nbins, hi, lo respectively.

If so, where do I fix the binning of this histo?

I think the problem starts with my roopdf. I am trying to fit the h1 histo.

TH1F *h1 = new TH1F("h1","e-e+ spectra;Mass;Events",200,0.,500.);
RooDataHist dh("dh","e-e+",x,Import(*h1)) ;

the dh histo has been drawn on the first panel and it seems to retain the count of the original histo (h1).

I am then fitting it with a pdf fullmodel

         RooRealVar x("x","mass",18.,500.);
         x.setRange("Range1",18.,75.) ;
         x.setRange("Range2",107.,495.) ;
         RooRealVar a1("a1","a1",-50,1000);
         RooRealVar a2("a2","a2",-50,400);
         RooRealVar a3("a3","a3",-50,400);
         RooRealVar a4("a4","a4",-500,1000);
         RooRealVar a5("a5","a5",-50,1000);
         RooRealVar a6("a6","a6",-50,1000);
         RooRealVar a7("a7","a7",-50,1000);
         RooRealVar a8("a8","a8",-50,100);
         RooRealVar a9("a9","a9",-50,100);
         RooRealVar aa("aa","aa",-50,100);  
         RooRealVar mass("mass","Central value of Gaussian",90.,80.,100.);
         RooRealVar sigma("sigma","Width of Gaussian",20,0,100);
         RooRealVar alpha("alpha","Alpha",40.,0.,100.);
         RooRealVar n("n","Order",6,0.,10.);

         RooCBShape crystalball("crystalball", "Signal Region", x,mass,sigma,alpha,n);
         RooRealVar b("b", "Number of background events", 0, 5500);
         RooRealVar s("s", "Number of signal events", 0, 500);

         RooBernstein bg_bern("bg_bern","background",x, RooArgList(a1,a2,a3,a4,a5,a6,a7,a8,a9));




          RooAddPdf fullModel("fullModel", "crystalball + bg_bern", RooArgList(crystalball, bg_bern),RooArgList(s, b));

         RooFitResult* r = fullModel.fitTo(dh,Save()) ;

But drawing this pdf on the second panel shows the low counts (because normalised?). This is where I am deriving the bg_bern pdf, and this must be where the low count issue arises.

TH1 *h2 = bg_bern.createHistogram("background;mass;count", x, Binning(200,0.,500.));

Still unable to populate the histogram properly

FYI, we’re in Christmas / New Year’s mode here - @StephanH might be able to help when he’s back, around mid-January.

I understand.
Merry Christmas and a happy new year to all of you! :smiley:

Hello,

you are right that you

  • have to match the binning of the histograms, so Binning(200, 0., 500.) is what you should do.
  • and you have to adjust the normalisation. Without doing anything, RooFit will always normalise the PDF. If you need custom scaling, you can either scale the histogram manually hxx->Add(hyy, -7000), or you can have it scaled when creating the histogram. Have a look at the documentation here:
    RooAbsPdf::createHistogram docs
    You will find the keyword Extended(), which instructs the function to normalise according to the event yield of the (extended) PDF. This only works if the PDF is an extended PDF, but that should be true in your case, as you are extending using the RooAddPdf.
    Let me know if Extended() doesn’t work.

One more thing
The RooAbsPdf, and therefore every PDF, has a function to retrieve the expected events:
RooAbsPdf::expectedEvents
This could be useful in your case to find the number of events you should scale the PDF to if you decide to do the scaling manually. You have to give it a RooArgSet(x) (all the variables of your PDF) such that it knows in which dimensions to compute the number of events.

Hello @StephanH,

Will that work with the bg_bern pdf too? I tried

TH1 *h2 = bg_bern.createHistogram("background;mass;count", x, Binning(200,0.,500.),Extended());

But that makes no difference

How do I create a RooArgSet out of the existing number of variables (mass,sigma,alpha,n, and a1 through a9) that I have?
I need to do

fullModel.expectedEvents(*rooargset goes here*)

No, because this one is not extended. You need to use the value of b to scale it, because b is exactly the number of events of the background component.

Actually, your RooArgSet would only be RooArgSet(x), because this is the only variable in your problem. All the others are parameters of your model. However, since you want the normalisation of the background component, you don’t need expectedEvents(), because this is exactly the value of b.

Is there a way I can convert a RooRealVar to a double?

I am trying h2 -> Scale( *double b goes here*) to scale the histo, as I believe the double value in the bracket is the new norm?

I can just hard code the number in, but is there a better way of going about it?

Reference : Different ways of normalizing histograms

Sure!
RooAbsReal::getVal()

Since I want the total (h1) - background (h2) histogram, I was doing h1 -> Add(h2, -1.0).

I want the updated h2 ( i.e. h2 counts + fullmodel.expectedEvents(RooArgSet(x)) which is b = 4017) subtracted from h1.

h2 -> Scale(4017.) gives me weird results. So I tried

h12 -> Add(h2, fullModel.expectedEvents(RooArgSet(x))) and that gives


Is this correct?

Is defining a new histogram recommended here? Something like a h12 to contain the scaled up h2. Then I can do h1 -> Add (h12, -1).
This looks messy, is there a better way of going about it?

What I meant here was:

h12 -> Add(h2, -1. * b.getVal())

This should scale the background component h2 to the number of background events that the fit extracts, and subtract it from the data.

Do you mean h12 -> Add(h2, 1*b.getVal()); ?

TH1 *h2 = bg_bern.createHistogram("background;mass;count", x, Binning(200,0.,500.));
cout << fullModel.expectedEvents(RooArgSet(x)) << endl;
TH1* h12 = new TH1I("h12", "scaled signal", 200, 0.0, 500.0);	
h12 -> Add(h2, 1.0*b.getVal());
		 
h1 -> Add(h12, -1.0);

I get


which looks okay except near 0

No, you want to subtract, don’t you? It looks like I got confused with h12, h1, h2 etc. It’s sufficient to subtract and scale directly. Is that not evident?

TH1 *h2 = bg_bern.createHistogram("background;mass;count", x, Binning(200,0.,500.));
h1 -> Add(h2, -1.0 * b.getVal());
1 Like