Looping over th2 to get the error in the histogram

iftikhar_safi · April 8, 2021, 10:37am

Hello everyone,
I have two th2 histograms, I want to divide one by another. However, the error bars(sqrt of entries) are too big if I do it this way. So I want to have the errors in each histogram calculated before, and then divide them too, the same way as I divide the histogram. The image shows the error after dividing one histogram by the other. The plot (th1) shown is the projection of each bin of the th2 histogram after division.
with parameters

I want to find a way to calculate the error bars for th2 before dividing it with another histogram, So that I do not have the errors bar too big?
Please, let me know if you need any more information.
Thank you for taking the time to read.

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

yus · April 8, 2021, 12:54pm

Hi Iftikhar,

so I guess you want to loop over all bins, to figure out two relative uncertainties on the bin content (one for the bin of the first histogram, another one for the same bin of the second histogram), add these two in quadrature, and assign this resulting relative uncertainty as an uncertainty on the ratio of the bin contents?

iftikhar_safi · April 8, 2021, 1:00pm

Hello Yus,
Thanks for replying,
I just want to know how can I find the errors(sqrt of entries) in th2 for a single th2 histogram. I just need the error is one 2-d histogram not the relative uncertainty. For th1 I did something like this

    for(int i=1; i<=h->GetNbinsX(); ++i) {
                        Double_t entry=h->GetBinContent(i);
                        Double_t error=TMath::Sqrt(entry);
                        h->SetBinContent(i,entry);
                        h->SetBinError(i, error);
                }

Is there any way I can do the same for a 2-d histogram?
Thanks

Regards
Iftikhar

yus · April 8, 2021, 1:04pm

Sure thing:

for(int i=1; i<=h->GetNbinsX(); ++i) {
     for(int j=1; j<=h->GetNbinsY(); ++j) {
         Double_t entry=h->GetBinContent(i, j);
         Double_t error=TMath::Sqrt(entry);
         h->SetBinContent(i, j, entry); // but what for? It's already set, isn't it?
         h->SetBinError(i, j, error);
    }
}

iftikhar_safi · April 8, 2021, 1:32pm

Hi Yus,
Thanks alot, I thought something like this will solve my problem.
However, I have the same results.

Let me explain more clearly if I may, I have a 2d histogram which contains number of counts in each day and there is another one which is livetime of the experiment in each day, both have the same number of days. I want to divide the 2d histogram which contains counts in each day by a 2d histogram which contains the livetime for each day.

Then after the division of these histogram I take projections of the new histogram, and then I have a peak in each projection, I integrate from a to b, and find the counts. The counts are then plotted in another histogram.

The problem is here (which is related to the error bars in the first th2 histograms which I used before the division), the error bars that I get are too big as you can see in the image I uploaded before…

So I want to find a way in which I have the error bar calculated for the 2d histogram which contains the counts in each day, and then I divide by the livetime in each day. The error bars are not important for the livetime in my case, just the error bars in counts in each day?

I hope its a little bit clear now.
Thank you very much.

Best,
Iftikhar

yus · April 8, 2021, 1:39pm

Why are the error bars too big? Were you expecting such large errors or was it due to a mistake somewhere? Of course, if you have large errors in histogram in the numerator, you’ll have large errors in the resulting ratio histogram. You should understand why the errors are too large in the initial histogram. If they should be like this, there is nothing you can do. If, on the other hand, it’s a mistake, then you should fix this when you were producing that histogram in the numerator, and not when you already have it.

iftikhar_safi · April 8, 2021, 1:46pm

Hello Yus,
Thanks
The error bars are big because I did not take into account the errors bars when I was dividing them.
Now I see that the errors bars are not right,so I want to take into account the error bars when I am dividing them.
I think if I can have the error bars on the th2 manually by taking the square root of the entries in each bin and then divide will solve the issue. I want to set the errors manually as shown below and then divide.

 Double_t error=TMath::Sqrt(entry);

Thank you very much.

Best,
Iftikhar

yus · April 8, 2021, 2:09pm

Ok, then you first reset the errors with

for(int i=1; i<=h->GetNbinsX(); ++i)
     for(int j=1; j<=h->GetNbinsY(); ++j)
         h->SetBinError(i, j, TMath::Sqrt(h->GetBinContent(i, j)));

and then you divide your histograms, right? Will this work for you?

iftikhar_safi · April 8, 2021, 2:18pm

Hello Yus,
Thanks for the code, However, It did not work for me. I don’t know if my approach is right.
I thought SetBinError will solve the issue. I do not know what can I do about it? How can I solve this issue any suggestions about this?

Thanks
Best,
Iftikhar

Wile_E_Coyote · April 8, 2021, 2:21pm

Provide a ROOT file, with both histograms that you want to divide, for inspection.

yus · April 8, 2021, 2:22pm

For starters, you could show us your histograms. I mean those you have in the very beginning, before you start manipulating with division and errors.

iftikhar_safi · April 8, 2021, 2:32pm

Hello,
Sure here are the three root files, The first one is counts in each day, the second one is livetime in each day and the third on is the division of counts by livetime.

events_days.root (71.4 KB)
Livetime_day.root (19.3 KB)
events_by_livetime.root (122.5 KB)

This is how they look like

I used the third one for the analysis i.e projections, counts and plotting them.
The analysis result is
with parameters

Thank you very much.

Best,
Iftikhar

yus · April 8, 2021, 2:38pm

But there are no histograms saved in these files (at least in those with days and lifetimes), are there? Can you upload the files with actual TH2 objects?

iftikhar_safi · April 8, 2021, 2:48pm

Hello,
Here is the script I am using for this.
normalized.C (14.8 KB)

The relevent part starts at line 166. Its a long script, sorry for that.
Please have a look at the link below
https://drive.google.com/drive/folders/1rsXXvgKtD2l3DtqYMRuE-cfHYlVGOgDr?usp=sharing
Thanks

Best,
Iftikhar

Wile_E_Coyote · April 8, 2021, 4:00pm

Comparing your “normalized.C” with your histograms … I don’t really understand what happens here …

Looking at the “hNePerRun_merged” from “events_days.root”, I can see huge errors, much bigger than “sqrt(bin_contents)”.
Try to print, e.g.:
hNePerRun_merged->GetBinContent(6, 6)
hNePerRun_merged->GetBinError(6, 6)

Looking at the “hLTmarged” from “Livetime_day.root”, I can see the same problem with huge errors.
Try to print, e.g.:
hLTmarged->GetBinContent(6, 6)
hLTmarged->GetBinError(6, 6)

Moreover, in the case of the “event counts”, I can understand that the errors should be the “sqrt” of them.
However, in the case of the “lifetimes”, I think the errors should be some percentage of them or some constant value (i.e., not the “sqrt”).

iftikhar_safi · April 8, 2021, 4:12pm

Hello Wile,

I do not understand what it shows, it give me two number with no relation to each other. Can you please give a brief explanation of what it means?
Thanks

Best,
Iftikhar

Wile_E_Coyote · April 8, 2021, 4:42pm

That’s what you need to explain (take this example bin and find out why its error is so big compared to its contents).

iftikhar_safi · April 12, 2021, 2:55pm

Hello,
I found why the error bars are wrong, I will explain briefly here.
I found the error bars by taking the square root of the counts and then dividing by the livetime to get the rate of the counts in a loop. Here the error(sqrt of counts) is the same as expected.

Now the problem is that in a new loop I take the projections of the 2d plot of rate of counts and I integrate from (a,b) in each bin, Then I plot the counts in a peak and my error bars goes up.

Is there any way to solve this issue, for example, If I can connect the error from the first loop to the counts that I get from the integral in the second loop?
Here is relevent script


   for (int ibin = 1; ibin <= h_counts->GetNbinsX(); ++ibin) { 
         for (int iybin = 1; iybin <= h_counts->GetNbinsY(); ++iybin) { 
              Double_t lvt = h_live->GetBinContent(ibin, iybin);
              Double_t entry = h_counts->GetBinContent(ibin, iybin);
              Double_t error = TMath::Sqrt(entry)/lvt;
              h_counts->SetBinContent(ibin, iybin, entry/lvt);
              h_counts->SetBinError(ibin, iybin, error);//
              cout<<"error: "<<h_counts->GetBinError(ibin, iybin)<<" sqrt:
             "<<TMath::Sqrt(h_counts->GetBinContent(ibin, iybin))<<endl;
                }
        }

and the next loop is below

  for (int64_t i = 1; i <=h_counts->GetNbinsX(); ++i) {
     TH1D *px_slices = h_counts->ProjectionY(Form("slices_%zu",i),i,i, "");
       px_slices->Add(background,-1);
       double counts_integral_peak1 =px_slices->Integral(xmin,xmax);
       h_rate->Fill(A,counts_integral_peak1);
       h_rate->Draw("e1");

How can I have the error from the first loop propagated to the second loop?
Thank you very much

Wile_E_Coyote · April 12, 2021, 4:38pm

TH1::IntegralAndError

iftikhar_safi · April 15, 2021, 9:43am

Hello,
Thank you very much Wile_E_Coyote.

Best
Iftikhar