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>
32 #include "AliUnicorHN.h"
36 //=============================================================================
37 AliUnicorHN::AliUnicorHN(Char_t *nam, Int_t ndim, TAxis **ax)
38 : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5),
41 // Constructor for building from scratch.
43 // Above, we just have managed to create a single one-dimensional histogram
44 // with number of bins equal to the product of the numbers of bins in all
45 // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5.
47 for (int i=0; i<fNdim; i++) ax[i]->SetName(Form("axis%d",i));
48 for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]);
50 // for speed reasons, number of bins of each axis is stored on fNbins
51 // and the running product of the latter on fMbins
52 // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1}
55 for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
56 for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
57 printf(" %d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX());
59 //=============================================================================
60 AliUnicorHN::AliUnicorHN(Char_t *filnam, Char_t *nam)
61 : TH1D(*((TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"))),
64 // Constructor for reading from file.
66 TFile *f = TFile::Open(filnam,"read");
68 TAxis *ax[fgkMaxNdim];
69 for (fNdim=0; fNdim<fgkMaxNdim; fNdim++) {
70 ax[fNdim] = (TAxis *) gDirectory->Get(Form("axis%d",fNdim));
71 if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]);
77 for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
78 for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
80 if (GetNbinsX()!=Albins(fNdim,ax)) {
81 printf("number of bins of histo %d differs from product of nbins of axes %d\n",
82 GetNbinsX(),Albins(fNdim,ax));
86 printf("%d-dimensional histogram %s with %d bins read from file %s\n",
87 fNdim,nam,GetNbinsX(),filnam);
89 //=============================================================================
90 Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax)
92 // Calculate product of nbins of ax[0]...ax[n-1]
93 // (= total number of bins of the multidimensional histogram to be made).
96 // while (n--) nbins *= ax[n]->GetNbins();
97 for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins();
100 //=============================================================================
101 Int_t AliUnicorHN::MulToOne(const Int_t * const k) const
103 // Calculate the 1-dim index n from n-dim indices k[fNdim].
104 // Valid k[i] should be between 0 and fNbins[i]-1.
105 // Valid result will be between 0 and GetNbinsX()-1.
106 // Return -1 if under- or overflow in any dimension.
109 for (int i=0; i<fNdim; i++) {
110 if (k[i]<0) return -1;
111 if (k[i]>=fNbins[i]) return -1;
116 //=============================================================================
117 Int_t AliUnicorHN::MulToOne(Double_t *x)
119 // Calculate the 1-dim index n from n-dim vector x, representing the
120 // abscissa of the n-dim histogram. The result will be between 0 and
124 for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1;
127 //=============================================================================
128 void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const
130 // Calculate the n-dim indices k[fNdim] from 1-dim index n.
131 // Valid n should be between 0 and GetNbinsX()-1.
132 // Valid results will be between 0 and fNbins[i]-1.
134 div_t idi; // integer division structure
135 for (int i=0; i<fNdim; i++) {
136 idi = div(n,fMbins[i]);
137 k[i] = idi.quot; // quotient
138 n = idi.rem; // remainder
141 //=============================================================================
142 Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w)
144 // Fill the histogram. The array xx holds the abscissa information, w is the
145 // weigth. The 1-dim histogram is filled using the standard Fill method
146 // (direct access to the arrays was tried and was not faster).
148 int nbin = MulToOne(xx);
149 if (nbin == -1) return 0;
150 return TH1D::Fill(nbin+1,w);
152 //=============================================================================
153 Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...)
155 // Fill the histogram. Arguments are passed as doubles rather than array.
157 Double_t xx[fgkMaxNdim] = {x0, x1};
159 for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t);
160 Double_t weigth = va_arg(ap,Double_t);
162 //for (int i=0; i<fgkMaxNdim; i++) printf("%8.2f ",xx[i]); printf("%6.2f \n",weigth);
163 return Fill(xx,weigth);
165 //=============================================================================
166 Int_t AliUnicorHN::Write() const
168 // Save the 1-dim histo and the axes in a subdirectory on file. This might
169 // not be the most elegant way but it is very simple and backward and forward
174 histo.SetName("histo"); // hadd bug fix; actually, does not cost memory, strange
176 TDirectory *dest = gDirectory->mkdir(GetName());
179 nbytes += histo.Write("histo");
180 for (int i=0; i<fNdim; i++) nbytes += fAxis[i].Write();
181 printf(" %s stored in %s\n",GetName(),dest->GetPath());
186 //=============================================================================
187 AliUnicorHN *AliUnicorHN::ProjectAlong(char *nam, Int_t dim, Int_t first, Int_t last)
189 // Reduce dimension dim by summing up its bins between first and last.
190 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
191 // Return the resulting fNdim-1 dimensional histogram.
193 if (dim<0 || dim>fNdim-1) return 0;
194 if (first<0) first = 1;
195 if (last<0) last = fNbins[dim];
197 // create new (reduced) histogram
199 TAxis *ax[fgkMaxNdim-1];
201 for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i);
202 char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
203 char *eame = Form("%s_error",name);
204 AliUnicorHN *his = new AliUnicorHN(name,fNdim-1,ax); // result histogram
205 AliUnicorHN *eis = new AliUnicorHN(eame,fNdim-1,ax); // temporary storage for errors
207 // sum up the content and errors squared
209 Int_t *k = new Int_t[fNdim]; // old hist multiindex
210 Int_t *m = new Int_t[fNdim-1]; // new hist multiindex
211 for (int i=0; i<GetNbinsX(); i++) {
213 if (k[dim]+1<first) continue;
214 if (k[dim]+1>last) continue;
216 for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j];
217 n = his->MulToOne(m);
218 his->AddBinContent(n+1,GetBinContent(i+1));
219 eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1));
222 // combine content and errors in one histogram
224 for (int i=0; i<his->GetNbinsX(); i++) {
225 his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1)));
228 his->SetLineColor(this->GetLineColor());
229 his->SetFillColor(this->GetFillColor());
230 his->SetMarkerColor(this->GetMarkerColor());
231 his->SetMarkerStyle(this->GetMarkerStyle());
238 //=============================================================================
239 TH1D *AliUnicorHN::ProjectOn(char *nam, Int_t dim, const Int_t * const first, const Int_t * const last) const
241 // Project on dimension dim. Use only bins between first[i] and last[i].
242 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
243 // first[i]=0 means from the first bin
244 // last[i]=0 means till the last bin
246 if (dim<0 || dim>fNdim-1) return 0;
248 // calculate the projection; lowest index 0
250 double *yy = new double[fNbins[dim]]; // value
251 double *ey = new double[fNbins[dim]]; // error
252 for (int i=0; i<fNbins[dim]; i++) yy[i]=0;
253 for (int i=0; i<fNbins[dim]; i++) ey[i]=0;
254 Int_t *k = new Int_t[fNdim];
255 for (int i=0; i<GetNbinsX(); i++) {
257 int isgood = 1; // bin within the range?
258 for (int j=0; j<fNdim; j++) {
259 if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
260 if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
263 yy[k[dim]]+=GetBinContent(i+1);
264 ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1);
268 // make the projection histogram
269 // use name nam if specified; otherwise generate one
272 char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
273 if (fAxis[dim].IsVariableBinSize())
274 his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
276 his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax());
277 his->SetXTitle(fAxis[dim].GetTitle());
278 // or his->GetXaxis()->ImportAttributes(ax);
280 his->SetLineColor(this->GetLineColor());
281 his->SetFillColor(this->GetFillColor());
282 his->SetMarkerColor(this->GetMarkerColor());
283 his->SetMarkerStyle(this->GetMarkerStyle());
284 for (int i=0; i<his->GetNbinsX(); i++) {
285 his->SetBinContent(i+1,yy[i]);
286 his->SetBinError(i+1,sqrt(ey[i]));
294 // if (name!=nam) delete [] name;
298 //=============================================================================
299 TH1D *AliUnicorHN::ProjectOn(char *nam, Int_t dim, const Double_t * const first, const Double_t * const last)
301 // Project on dimension dim. Use only bins between first[i] and last[i].
303 if (dim<0 || dim>fNdim-1) return 0;
304 Int_t kfirst[fgkMaxNdim];
305 Int_t klast[fgkMaxNdim];
306 for (int i=0; i<fNdim; i++) {
307 kfirst[i] = fAxis[i].FindBin(first[i]);
308 klast[i] = fAxis[i].FindBin(last[i]);
310 return ProjectOn(nam,dim,kfirst,klast);
312 //=============================================================================
313 TH2D *AliUnicorHN::ProjectOn(char *nam, Int_t dim0, Int_t dim1, const Int_t * const first, const Int_t * const last) const
315 // Project on dim1 vs dim0. Use only bins between first[i] and last[i].
316 // Use root convention: bin=1 is the first bin, bin=nbins is the last.
317 // first[i]=0 means from the first bin
318 // last[i]=0 means till the last bin
320 if (dim0<0 || dim0>fNdim-1) return 0;
321 if (dim1<0 || dim1>fNdim-1) return 0;
323 // calculate the projection
325 double **yy = new double*[fNbins[dim0]]; // value
326 double **ey = new double*[fNbins[dim0]]; // error
327 for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]];
328 for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]];
329 for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0;
330 for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0;
331 Int_t *k = new Int_t[fNdim];
332 for (int i=0; i<GetNbinsX(); i++) {
334 int isgood = 1; // bin within the range?
335 for (int j=0; j<fNdim; j++) {
336 if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
337 if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
340 yy[k[dim0]][k[dim1]]+=GetBinContent(i+1);
341 ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1);
345 // make the projection histogram
346 // use name nam if specified; otherwise generate one
349 char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0);
350 if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
351 his = new TH2D(name,name,
352 fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
353 fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
354 else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
355 his = new TH2D(name,name,
356 fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
357 fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
358 else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
359 his = new TH2D(name,name,
360 fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
361 fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
362 else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
363 his = new TH2D(name,name,
364 fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
365 fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
366 his->SetXTitle(fAxis[dim0].GetTitle());
367 his->SetYTitle(fAxis[dim1].GetTitle());
368 // or his->GetXaxis()->ImportAttributes(ax0);
369 // or his->GetYaxis()->ImportAttributes(ax1);
371 his->SetLineColor(this->GetLineColor());
372 his->SetFillColor(this->GetFillColor());
373 his->SetMarkerColor(this->GetMarkerColor());
374 his->SetMarkerStyle(this->GetMarkerStyle());
375 for (int i=0; i<his->GetNbinsX(); i++)
376 for (int j=0; j<his->GetNbinsY(); j++) {
377 his->SetBinContent(i+1,j+1,yy[i][j]);
378 his->SetBinError(i+1,j+1,sqrt(ey[i][j]));
383 for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i];
384 for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i];
388 // if (name!=nam) delete [] name;
392 //=============================================================================