I have a bunch of (say, 100) 2D histograms of equal dimensions. Now I would like to reduce their information to two histograms: one where each bin represents the average of the same bin from all the 100 histograms, and a second one, where each bin represents the standard deviation of that bin from all histograms. Point one is easy with TH1::Add and then division by the number of histograms, but what would be the most efficient way to get to point two (standard deviation for each bin)? Since my histograms/arrays are large (order 1000x1000), an efficient computation is very important (but I will also be grateful for simple, inefficient solutions to begin with).
OK, I had a look at the documentation. While TProfile2D is likely what I need, it is not clear to me how to construct the TProfile2D from a list of histograms (other than filling the entries “by hand”). Do I have to loop over the bin contents or is there a more efficient way?
In the following simple example, what is the missing code to go from the TList of TH2D to the TProfile2D, such that the TProfile2D will contain bin-by-bin average and RMS of the 100 histograms? I cannot use TProfile2D::Merge()…
{
Int_t nBinX = 5;
Int_t nBinY = 5;
Double_t xMin = 0.;
Double_t xMax = 100.;
Double_t yMin = 0.;
Double_t yMax = 100.;
Int_t nHist = 100;
// fill 100 histograms
TH2D* hist[100];
for (Int_t i = 0; i < nHist; i++) {
TString name = TString("hist");
name += i;
hist[i] = new TH2D(name.Data(),"hist", nBinX, xMin, xMax, nBinY, yMin, yMax);
hist[i]->FillRandom("expo",10000);
}
// put them into a TList:
TList *list = new TList;
for (Int_t i = 0;i <nHist; i++) {
list->Add(hist[i]);
}
// combine them to a TProfile2D?
TCanvas *c1 = new TCanvas("c1","Profile histogram example",200,10,700,500);
hprof2d = new TProfile2D();
??? hprof2d->IsThereSuchAMethod(list); ???
hprof2d->Draw("box");
}
It seems like there is no simple way to extract the TProfile2D from a bunch of histograms. For reference, I include here an example demonstrating the operation with loops over the histogram entries. Would it be worthwhile to encapsulate this into a method of TProfile2D or am I the only one having this problem?
// Example for how to merge a bunch of 2D histograms into one which contains bin-by-bin
// average and RMS of the individual histogram bins.
{
Int_t nBinX = 5;
Int_t nBinY = 5;
Double_t xMin = 0.;
Double_t xMax = 1.;
Double_t yMin = 0.;
Double_t yMax = 1.;
Int_t nHist = 100;
// create and fill 100 histograms
TH2D* hist[100];
for (Int_t i = 0; i < nHist; i++) {
TString name = TString("hist");
name += i;
hist[i] = new TH2D(name.Data(),"hist", nBinX, xMin, xMax, nBinY, yMin, yMax);
hist[i]->FillRandom("expo",100000);
}
// combine them to a profile histogram containing average and rms, bin by bin:
hprof2d = new TProfile2D("hprof2d","hprof2d", nBinX, xMin, xMax, nBinY, yMin, yMax,"s");
for (Int_t i = 0; i < nHist; i++) {
for (Int_t ix = 1; ix <= nBinX; ix++) {
for (Int_t iy = 1; iy <= nBinY; iy++) {
Double_t x = hist[i]->GetXaxis()->GetBinCenter(ix);
Double_t y = hist[i]->GetYaxis()->GetBinCenter(iy);
hprof2d->Fill(x,y,hist[i]->GetBinContent(ix,iy),1);
}
}
}
hprof2d->Draw("box");
// check content of bin (2,3):
TGraph* gr = new TGraph(100);
for (Int_t i = 0; i < nHist; i++) {
cout << "histogram nr. " << i << " , content of bin (2,3) is " << hist[i]->GetBinContent(2,3) << endl;
gr->SetPoint(i,hist[i]->GetBinContent(2,3),1.);
}
cout << endl;
cout << " profile histo: content of bin (2,3) is " << hprof2d->GetBinContent(2,3) << endl;
cout << " profile histo: error of bin (2,3) is " << hprof2d->GetBinError(2,3) << endl;
cout << " mean from graph: " << gr->GetMean(1) << endl;
cout << " rms from graph: " << gr->GetRMS(1) << endl;
}
The histograms are given. They are images obtained from a camera. The idea then is to take 100 images and compress this information into one histogram that has average and rms deviation (for background subtraction).
I think that it is difficult to justify the addition of a new function in TProfile2D just to replace a few lines of code in your very particular example. By the way, I suggest a simplified and faster loop with:
// combine them to a profile histogram containing average and rms, bin by bin:
hprof2d = new TProfile2D("hprof2d","hprof2d", nBinX, xMin, xMax, nBinY, yMin, yMax,"s");
for (Int_t ix = 1; ix <= nBinX; ix++) {
Double_t x = hprof2d->GetXaxis()->GetBinCenter(ix);
for (Int_t iy = 1; iy <= nBinY; iy++) {
Double_t y = hprof2d->GetYaxis()->GetBinCenter(iy);
for (Int_t i = 0; i < nHist; i++) {
hprof2d->Fill(x,y,hist[i]->GetBinContent(ix,iy),1);
}
}
}