Home | News | Documentation | Download

Covariance of two histograms

Hello,

I have two histograms.
Like,

TH1D* a=...; 
TH1D* b=...;

I want to get the covariance bewtween a and b. Is there a simple method to get
the covariance without recreate a TH2D( );
Thanks.

Didn’t you already posted the same question (the one you just withdrawn)? Didn’t @moneta already answered the question?

@bellenot, Hi bellenot,
Yes, @moneta has helped me resolve how to propagate a uncertainty between two variables.
Actually, to get the uncertainty, a covariance between two variables is necessary.
I get a similar question like below. However, I need to recreate a new TH2. I wonder whether there is
a more simple method to get it.

Thanks a lot.
Liuyao

Hi

I though you need to create a covariance between each single bin a(i) bs b(i) for used for adding histograms as you asked here:

The function TH2::GetCovariance return the global correlation between all values of a and b.

You would need to compute for each single bin (at filling time) something shown here

Lorenzo

@moneta,
Thanks, could you please make a simple example for me?.
Actually, I am a little confused to create a covariance between each single bin a[i] , b[i].

Hi,
I can provide you an example, but it is maybe better you clarify exactly what is your problem, what are your data and what do you want to do, and then I can provide you an example which will be useful for you

Lorenzo

Hello @moneta,
I am sorry to reply to you so late, due to something other occupied.
Here, I list my problem exactly.
I want to calculate a sum of mean a and mean b vs mass in all of event, namely

(<<a>> - <<b>>) vs mass

so I created two TH2D:

a->Fill(mass, <a>*weight1),
b->Fill(mass, <b>*weight2),   

Next part is the calculation process,

TH2D* a = (TH2D*)file1->Get(" ");
TH2D* b = (TH2D*)file1->Get(" ");
TH2D *acopy = (TH2D*) a->Clone();

TProfile*  tp_a = (TProfile*)acopy->ProfileX("tp_a");
TProfile*  tp_b = (TPrifile*)b->ProfileX("tp_b");

TH1D* th_a=(TH1D*) tp_a->ProjectionX("");
TH1D* th_b=(TH1D*) tp_b->ProjectionX("");

th_a->Add(th_b, -1); 

Actually, <<a>> and <<b>>, they are correlated, not independent.
So I have to consider the error of sum.
I know they have Reliability weights.
Refer to the below formula. <<a>>, <<b>> representing to <<4>>,<<2>>

Firstly, I have a problems.

  1. To get ∑(w^2), here, I consider w is the weight in each TProfile bin. So I want to use GetBinSumw2() , But the GetBinSumw2 return an array, I want to get a value in each bin.
  2. How to calcualte the Cov(<2>,<4>).

I am appreciate if you can give me a example.

Hi,
What are weight1 and weight2 ? If these are the weights you should fill the 2D histogram as this:

hist_a->Fill(mass, var_a, weight1);
hist_b->Fill(mass,  var_b, weight2); 

Where I call var_a the sample values of the variable a that you have.

Do you need to compute the variance of QC, s2(QC) for each bin mass value ?

Lorenzo

Hi Moneta,
Yes, you are right. Both a and b are TH2D, and they are filled Fill(mass, <2>/weight, weight )
//here, weight represents the multiplicity.
I uplated a sample include two TH2D <2> vs mass distribution (hist_Q2Q2_D0plusD0barHF_Re_cent0to10_pT1.0to10.0_y-1.0to1.0) and <4> vs mass distribution (hist_Q2Q2Q2Q2_D0plusD0barHFHFHF_Re_cent0to10_pT1.0to10.0_y-1.0to1.0) .

Yes, my purpose to calcuate the ‘QC, s2(QC)’, I have calculate the QC. So I just confused to the
calculation of s2(QC).

Thanks.

And sample.root (33.8 KB)

@moneta,
Hi Moneta,

For the calculation process of s2QC{4},
Through the analysis I posted before, we know <<a>>, <<b>> representing the content of TH1D ( x representing mass). I can calculate first and second part of the formula:


use below codes:

Double_t stats[7];
Double_t sumw_a; 
Double_t sumw2_a; 
Double_t sumw_b; 
Double_t sumw2_b; 
for(Int_t i=1; i<th_a->GetNbinsX()+1; i++){
    th_a->GetXaxis()->SetRange(i, i);
    th_a->GetStats(stats);
    sumw_a  = stats[0]; //sumw in each bin; 
    sumw2_a = stats[1]; //sumw2 in each bin; 
    16* pow(th_a->GetBinContent(i), 2) * sumw2_a /(sumw_a*sumw_a) * pow(th_a->GetBinerror(i),2); //first part.  
    sumw2_b/(sumw_b*sumw_b) * pow(th_b->GetBinError(i),2); //second part.  

How for the third part of the formula, I don’t know how to calculate ∑w1∑w2 and cov(<2>,<4>).
I am appreciated if you can help me.
Best wishes.
Liuyao