1 /*************************************************************************
2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
18 //=============================================================================
19 // multidimensional histogram
20 // Physically, the data are kept in one single one-dimensional histogram.
21 // Projecting on 1, 2, and n-1 dimensions is implemented.
22 // The histogram can be saved on root file as the one-dimensional data
23 // histogram and the axes, thus eternal forward compatibility is ensured.
24 //=============================================================================
29 #include <TDirectory.h>
33 #include <TObjArray.h>
34 #include <TObjString.h>
36 #include "AliUnicorHN.h"
40 //=============================================================================
41 AliUnicorHN::AliUnicorHN(const char *nam, Int_t ndim, TAxis **ax)
42 : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5),
47 // Above, we just have managed to create a single one-dimensional histogram
48 // with number of bins equal to the product of the numbers of bins in all
49 // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5.
51 for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]);
52 for (int i=0; i<fNdim; i++) fAxis[i].SetName(Form("axis%d",i));
54 // for speed reasons, number of bins of each axis is stored on fNbins
55 // and the running product of the latter on fMbins
56 // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1}
58 for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
59 for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
60 for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
61 printf(" %d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX());
63 //=============================================================================
64 AliUnicorHN* AliUnicorHN::Retrieve(const char *filnam, const char *nam) {
66 // retrieve a multidimensional histogram from file
68 TFile *f = TFile::Open(filnam,"read");
70 if (!f->cd(nam)) {f->Close(); return 0;}
71 TH1D *hist = (TH1D*) gDirectory->Get("histo");
72 if (!hist) {printf("No histogram histo on file %s in directory %s\n",filnam,nam); return 0;}
73 hist->SetDirectory(gROOT);
74 TAxis *ax[fgkMaxNdim]={0};
76 while ((ax[n] = (TAxis *) gDirectory->Get(Form("axis%d",n)))) n++;
79 if (hist->GetNbinsX()!=Albins(n,ax)) {
80 printf("number of bins of histo %d differs from product of nbins of axes %d\n",
81 hist->GetNbinsX(),Albins(n,ax));
85 // derive the name from the path (ccc from aaa/bbb/ccc)
87 TObjArray *ar = path.Tokenize("/");
88 TObjString *os = (TObjString*) ar->Last();
89 const char *lastnam = os->GetString().Data();
91 AliUnicorHN *hi = new AliUnicorHN(lastnam,n,ax);
92 // *((TH1D*) hi) = *hist;
93 hist->Copy(*((TH1D*)hi));
98 //=============================================================================
99 Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax) {
101 // Calculate product of nbins of ax[0]...ax[n-1]
102 // (= total number of bins of the multidimensional histogram to be made).
105 // while (n--) nbins *= ax[n]->GetNbins();
106 for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins();
109 //=============================================================================
110 Int_t AliUnicorHN::MulToOne(const Int_t * const k) const {
112 // Calculate the 1-dim index n from n-dim indices k[fNdim].
113 // Valid k[i] should be between 0 and fNbins[i]-1.
114 // Valid result will be between 0 and GetNbinsX()-1.
115 // Return -1 if under- or overflow in any dimension.
118 for (int i=0; i<fNdim; i++) {
119 if (k[i]<0) return -1;
120 if (k[i]>=fNbins[i]) return -1;
125 //=============================================================================
126 Int_t AliUnicorHN::MulToOne(Double_t *x) {
128 // Calculate the 1-dim index n from n-dim vector x, representing the
129 // abscissa of the n-dim histogram. The result will be between 0 and
132 Int_t k[fgkMaxNdim] = {0};
133 for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindFixBin(x[i])-1;
136 //=============================================================================
137 void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const {
139 // Calculate the n-dim indices k[fNdim] from 1-dim index n.
140 // Valid n should be between 0 and GetNbinsX()-1.
141 // Valid results will be between 0 and fNbins[i]-1.
143 div_t idi; // integer division structure
144 for (int i=0; i<fNdim; i++) {
145 idi = div(n,fMbins[i]);
146 k[i] = idi.quot; // quotient
147 n = idi.rem; // remainder
150 //=============================================================================
151 Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w) {
153 // Fill the histogram. The array xx holds the abscissa information, w is the
154 // weigth. The 1-dim histogram is filled using the standard Fill method
155 // (direct access to the arrays was tried and was not faster).
157 int nbin = MulToOne(xx);
158 if (nbin == -1) return 0;
159 return TH1D::Fill(nbin+1,w);
161 //=============================================================================
162 Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...) {
164 // Fill the histogram. Arguments are passed as doubles rather than array.
167 Double_t xx[fgkMaxNdim] = {x0, x1};
169 for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t);
170 Double_t weigth = va_arg(ap,Double_t);
172 return Fill(xx,weigth);
174 //=============================================================================
175 Int_t AliUnicorHN::Save() const {
177 // Save the 1-dim histo and the axes in a subdirectory on file. This might
178 // not be the most elegant way but it is very simple and backward and forward
183 histo.SetName("histo"); // hadd bug fix; actually, does not cost memory, strange
185 TDirectory *dest = gDirectory->mkdir(GetName());
188 nbytes += histo.Write("histo");
189 for (int i=0; i<fNdim; i++) nbytes += fAxis[i].Write();
190 printf(" %s stored in %s\n",GetName(),dest->GetPath());
195 //=============================================================================
196 AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first, Int_t last) {
198 // Reduce dimension dim by summing up its bins between first and last.
199 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
200 // last=0 means till the last bin
201 // Return the resulting fNdim-1 dimensional histogram.
203 if (dim<0 || dim>fNdim-1) return 0;
204 if (last<=0) last = fNbins[dim];
206 // create new (reduced) histogram
208 TAxis *ax[fgkMaxNdim-1] = {0};
210 for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i);
211 const char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
212 const char *eame = Form("%s_error",name);
213 AliUnicorHN *his = new AliUnicorHN(name,fNdim-1,ax); // result histogram
214 AliUnicorHN *eis = new AliUnicorHN(eame,fNdim-1,ax); // temporary storage for errors
216 // sum up the content and errors squared
218 int k[fgkMaxNdim] = {0}; // old hist multiindex
219 int m[fgkMaxNdim] = {0}; // new hist multiindex
220 for (int i=0; i<GetNbinsX(); i++) {
222 if (k[dim]+1<first) continue;
223 if (k[dim]+1>last) continue;
225 for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j];
226 n = his->MulToOne(m);
227 his->AddBinContent(n+1,GetBinContent(i+1));
228 eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1));
231 // combine content and errors in one histogram
233 for (int i=0; i<his->GetNbinsX(); i++) {
234 his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1)));
237 his->SetLineColor(this->GetLineColor());
238 his->SetFillColor(this->GetFillColor());
239 his->SetMarkerColor(this->GetMarkerColor());
240 his->SetMarkerStyle(this->GetMarkerStyle());
247 //=============================================================================
248 TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first,
249 const Int_t * const last) const {
251 // Project on dimension dim. Use only bins between first[i] and last[i].
252 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
253 // first[i]=0 means from the first bin
254 // last[i]=0 means till the last bin
256 if (dim<0 || dim>fNdim-1) return 0;
258 // calculate the projection; lowest index 0
260 double *yy = new double[fNbins[dim]]; // value
261 double *ey = new double[fNbins[dim]]; // error
262 for (int i=0; i<fNbins[dim]; i++) yy[i]=0;
263 for (int i=0; i<fNbins[dim]; i++) ey[i]=0;
264 Int_t *k = new Int_t[fNdim];
265 for (int i=0; i<GetNbinsX(); i++) {
267 int isgood = 1; // bin within the range?
268 for (int j=0; j<fNdim; j++) {
269 if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
270 if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
273 yy[k[dim]]+=GetBinContent(i+1);
274 ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1);
278 // make the projection histogram
279 // use name nam if specified; otherwise generate one
282 const char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
283 if (fAxis[dim].IsVariableBinSize())
284 his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
286 his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax());
287 his->SetXTitle(fAxis[dim].GetTitle());
288 // or his->GetXaxis()->ImportAttributes(ax);
290 his->SetLineColor(this->GetLineColor());
291 his->SetFillColor(this->GetFillColor());
292 his->SetMarkerColor(this->GetMarkerColor());
293 his->SetMarkerStyle(this->GetMarkerStyle());
294 for (int i=0; i<his->GetNbinsX(); i++) {
295 his->SetBinContent(i+1,yy[i]);
296 his->SetBinError(i+1,sqrt(ey[i]));
304 // if (name!=nam) delete [] name;
308 //=============================================================================
309 TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first,
310 const Double_t * const last) const {
312 // Project on dimension dim. Use only bins between first[i] and last[i].
314 if (dim<0 || dim>fNdim-1) return 0;
315 Int_t kfirst[fgkMaxNdim] = {0};
316 Int_t klast[fgkMaxNdim] = {0};
317 for (int i=0; i<fNdim; i++) {
318 kfirst[i] = fAxis[i].FindFixBin(first[i]);
319 klast[i] = fAxis[i].FindFixBin(last[i]);
321 return ProjectOn(nam,dim,kfirst,klast);
323 //=============================================================================
324 TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t *
325 const first, const Int_t * const last) const {
327 // Project on dim1 vs dim0. Use only bins between first[i] and last[i].
328 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
329 // first[i]=0 means from the first bin
330 // last[i]=0 means till the last bin
332 if (dim0<0 || dim0>fNdim-1) return 0;
333 if (dim1<0 || dim1>fNdim-1) return 0;
335 // calculate the projection
337 double **yy = new double*[fNbins[dim0]]; // value
338 double **ey = new double*[fNbins[dim0]]; // error
339 for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]];
340 for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]];
341 for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0;
342 for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0;
343 Int_t *k = new Int_t[fNdim];
344 for (int i=0; i<GetNbinsX(); i++) {
346 int isgood = 1; // bin within the range?
347 for (int j=0; j<fNdim; j++) {
348 if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
349 if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
352 yy[k[dim0]][k[dim1]]+=GetBinContent(i+1);
353 ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1);
357 // make the projection histogram
358 // use name nam if specified; otherwise generate one
361 const char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0);
362 if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
363 his = new TH2D(name,name,
364 fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
365 fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
366 else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
367 his = new TH2D(name,name,
368 fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
369 fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
370 else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
371 his = new TH2D(name,name,
372 fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
373 fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
374 else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
375 his = new TH2D(name,name,
376 fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
377 fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
378 his->SetXTitle(fAxis[dim0].GetTitle());
379 his->SetYTitle(fAxis[dim1].GetTitle());
380 // or his->GetXaxis()->ImportAttributes(ax0);
381 // or his->GetYaxis()->ImportAttributes(ax1);
383 his->SetLineColor(this->GetLineColor());
384 his->SetFillColor(this->GetFillColor());
385 his->SetMarkerColor(this->GetMarkerColor());
386 his->SetMarkerStyle(this->GetMarkerStyle());
387 for (int i=0; i<his->GetNbinsX(); i++)
388 for (int j=0; j<his->GetNbinsY(); j++) {
389 his->SetBinContent(i+1,j+1,yy[i][j]);
390 his->SetBinError(i+1,j+1,sqrt(ey[i][j]));
395 for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i];
396 for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i];
400 // if (name!=nam) delete [] name;
404 //=============================================================================
405 TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Double_t *
406 const first, const Double_t * const last) const {
408 // Project on dim1 vs dim0. Use only bins between first[i] and last[i].
410 if (dim0<0 || dim0>fNdim-1) return 0;
411 if (dim1<0 || dim1>fNdim-1) return 0;
412 Int_t kfirst[fgkMaxNdim] = {0};
413 Int_t klast[fgkMaxNdim] = {0};
414 for (int i=0; i<fNdim; i++) {
415 kfirst[i] = fAxis[i].FindFixBin(first[i]);
416 klast[i] = fAxis[i].FindFixBin(last[i]);
418 return ProjectOn(nam,dim0,dim1,kfirst,klast);
420 //=============================================================================