Subtraction of two roodatahists

Debajyoti_Sengupta · December 16, 2018, 10:57am

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).

StephanH · December 17, 2018, 7:37pm

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.

Debajyoti_Sengupta · December 21, 2018, 1:47pm

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.?

Wile_E_Coyote · December 21, 2018, 6:15pm

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();

Debajyoti_Sengupta · December 21, 2018, 6:34pm

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)

Wile_E_Coyote · December 21, 2018, 6:43pm

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).

Debajyoti_Sengupta · December 22, 2018, 3:13am

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?

Debajyoti_Sengupta · December 22, 2018, 4:22am

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.));

Debajyoti_Sengupta · December 25, 2018, 12:37pm

Still unable to populate the histogram properly

Axel · December 31, 2018, 11:20am

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

Debajyoti_Sengupta · December 31, 2018, 11:23am

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

StephanH · January 8, 2019, 8:37am

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.

Debajyoti_Sengupta · January 8, 2019, 9:12am

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*)

StephanH · January 8, 2019, 1:02pm

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.

Debajyoti_Sengupta:

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*)

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.

Debajyoti_Sengupta · January 8, 2019, 2:27pm

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

StephanH · January 8, 2019, 3:58pm

Sure!
RooAbsReal::getVal()

Debajyoti_Sengupta · January 8, 2019, 4:05pm

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?

StephanH · January 8, 2019, 4:20pm

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.

Debajyoti_Sengupta · January 8, 2019, 4:26pm

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

StephanH · January 9, 2019, 9:20am

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());