]>
Commit | Line | Data |
---|---|---|
621688e4 | 1 | /************************************************************************* |
2 | * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
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 | **************************************************************************/ | |
15 | ||
16 | // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007 | |
17 | ||
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 | //============================================================================= | |
25 | ||
26 | #include <cmath> | |
27 | #include <stdlib.h> | |
28 | #include <TFile.h> | |
29 | #include <TDirectory.h> | |
28eee19b | 30 | #include <TROOT.h> |
621688e4 | 31 | #include <TAxis.h> |
32 | #include <TH2.h> | |
28eee19b | 33 | #include <TObjArray.h> |
34 | #include <TObjString.h> | |
35 | #include <TString.h> | |
621688e4 | 36 | #include "AliUnicorHN.h" |
37 | ||
38 | ClassImp(AliUnicorHN) | |
39 | ||
40 | //============================================================================= | |
61e4657c | 41 | AliUnicorHN::AliUnicorHN(const char *nam, Int_t ndim, TAxis **ax) |
621688e4 | 42 | : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5), |
28eee19b | 43 | fNdim(ndim) { |
44 | ||
45 | // constructor | |
621688e4 | 46 | |
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. | |
50 | ||
51 | for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]); | |
28eee19b | 52 | for (int i=0; i<fNdim; i++) fAxis[i].SetName(Form("axis%d",i)); |
53 | ||
621688e4 | 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} | |
28eee19b | 57 | |
61e4657c | 58 | for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins(); |
59 | for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1; | |
621688e4 | 60 | for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i]; |
c6fc7f72 | 61 | printf(" %d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX()); |
621688e4 | 62 | } |
63 | //============================================================================= | |
28eee19b | 64 | AliUnicorHN* AliUnicorHN::Retrieve(const char *filnam, const char *nam) { |
621688e4 | 65 | |
28eee19b | 66 | // retrieve a multidimensional histogram from file |
621688e4 | 67 | |
28eee19b | 68 | TFile *f = TFile::Open(filnam,"read"); |
69 | if (!f) return 0; | |
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}; | |
75 | int n=0; | |
76 | while ((ax[n] = (TAxis *) gDirectory->Get(Form("axis%d",n)))) n++; | |
77 | f->Close(); | |
78 | ||
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)); | |
82 | return 0; | |
621688e4 | 83 | } |
28eee19b | 84 | |
85 | // derive the name from the path (ccc from aaa/bbb/ccc) | |
86 | TString path=nam; | |
87 | TObjArray *ar = path.Tokenize("/"); | |
88 | TObjString *os = (TObjString*) ar->Last(); | |
89 | const char *lastnam = os->GetString().Data(); | |
90 | ||
91 | AliUnicorHN *hi = new AliUnicorHN(lastnam,n,ax); | |
92 | // *((TH1D*) hi) = *hist; | |
93 | hist->Copy(*((TH1D*)hi)); | |
94 | hi->SetName(lastnam); | |
95 | delete hist; | |
09d5920f | 96 | delete ar; |
28eee19b | 97 | return hi; |
621688e4 | 98 | } |
99 | //============================================================================= | |
28eee19b | 100 | Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax) { |
101 | ||
621688e4 | 102 | // Calculate product of nbins of ax[0]...ax[n-1] |
103 | // (= total number of bins of the multidimensional histogram to be made). | |
104 | ||
105 | Int_t nbins = 1; | |
106 | // while (n--) nbins *= ax[n]->GetNbins(); | |
107 | for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins(); | |
108 | return nbins; | |
109 | } | |
110 | //============================================================================= | |
28eee19b | 111 | Int_t AliUnicorHN::MulToOne(const Int_t * const k) const { |
112 | ||
621688e4 | 113 | // Calculate the 1-dim index n from n-dim indices k[fNdim]. |
114 | // Valid k[i] should be between 0 and fNbins[i]-1. | |
115 | // Valid result will be between 0 and GetNbinsX()-1. | |
116 | // Return -1 if under- or overflow in any dimension. | |
117 | ||
118 | Int_t n = 0; | |
119 | for (int i=0; i<fNdim; i++) { | |
120 | if (k[i]<0) return -1; | |
121 | if (k[i]>=fNbins[i]) return -1; | |
122 | n += fMbins[i]*k[i]; | |
123 | } | |
124 | return n; | |
125 | } | |
126 | //============================================================================= | |
28eee19b | 127 | Int_t AliUnicorHN::MulToOne(Double_t *x) { |
128 | ||
621688e4 | 129 | // Calculate the 1-dim index n from n-dim vector x, representing the |
130 | // abscissa of the n-dim histogram. The result will be between 0 and | |
131 | // GetNbinsX()-1. | |
132 | ||
61e4657c | 133 | Int_t k[fgkMaxNdim] = {0}; |
28eee19b | 134 | for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindFixBin(x[i])-1; |
621688e4 | 135 | return MulToOne(k); |
136 | } | |
137 | //============================================================================= | |
28eee19b | 138 | void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const { |
139 | ||
621688e4 | 140 | // Calculate the n-dim indices k[fNdim] from 1-dim index n. |
141 | // Valid n should be between 0 and GetNbinsX()-1. | |
142 | // Valid results will be between 0 and fNbins[i]-1. | |
143 | ||
144 | div_t idi; // integer division structure | |
145 | for (int i=0; i<fNdim; i++) { | |
146 | idi = div(n,fMbins[i]); | |
147 | k[i] = idi.quot; // quotient | |
148 | n = idi.rem; // remainder | |
149 | } | |
150 | } | |
151 | //============================================================================= | |
28eee19b | 152 | Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w) { |
153 | ||
621688e4 | 154 | // Fill the histogram. The array xx holds the abscissa information, w is the |
155 | // weigth. The 1-dim histogram is filled using the standard Fill method | |
156 | // (direct access to the arrays was tried and was not faster). | |
157 | ||
158 | int nbin = MulToOne(xx); | |
159 | if (nbin == -1) return 0; | |
160 | return TH1D::Fill(nbin+1,w); | |
161 | } | |
162 | //============================================================================= | |
28eee19b | 163 | Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...) { |
164 | ||
621688e4 | 165 | // Fill the histogram. Arguments are passed as doubles rather than array. |
28eee19b | 166 | |
621688e4 | 167 | va_list ap; |
168 | Double_t xx[fgkMaxNdim] = {x0, x1}; | |
169 | va_start(ap,x1); | |
170 | for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t); | |
171 | Double_t weigth = va_arg(ap,Double_t); | |
172 | va_end(ap); | |
173 | return Fill(xx,weigth); | |
174 | } | |
175 | //============================================================================= | |
28eee19b | 176 | Int_t AliUnicorHN::Save() const { |
177 | ||
621688e4 | 178 | // Save the 1-dim histo and the axes in a subdirectory on file. This might |
179 | // not be the most elegant way but it is very simple and backward and forward | |
180 | // compatible. | |
181 | ||
182 | Int_t nbytes = 0; | |
183 | TH1D histo(*this); | |
184 | histo.SetName("histo"); // hadd bug fix; actually, does not cost memory, strange | |
185 | ||
186 | TDirectory *dest = gDirectory->mkdir(GetName()); | |
187 | if (dest) { | |
188 | dest->cd(); | |
189 | nbytes += histo.Write("histo"); | |
c6fc7f72 | 190 | for (int i=0; i<fNdim; i++) nbytes += fAxis[i].Write(); |
621688e4 | 191 | printf(" %s stored in %s\n",GetName(),dest->GetPath()); |
192 | dest->cd(".."); | |
193 | } | |
194 | return nbytes; | |
195 | } | |
196 | //============================================================================= | |
28eee19b | 197 | AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first, Int_t last) { |
198 | ||
621688e4 | 199 | // Reduce dimension dim by summing up its bins between first and last. |
200 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
76d78859 | 201 | // last=0 means till the last bin |
621688e4 | 202 | // Return the resulting fNdim-1 dimensional histogram. |
203 | ||
204 | if (dim<0 || dim>fNdim-1) return 0; | |
76d78859 | 205 | if (last<=0) last = fNbins[dim]; |
621688e4 | 206 | |
207 | // create new (reduced) histogram | |
208 | ||
61e4657c | 209 | TAxis *ax[fgkMaxNdim-1] = {0}; |
621688e4 | 210 | int n=0; |
211 | for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i); | |
61e4657c | 212 | const char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim); |
213 | const char *eame = Form("%s_error",name); | |
621688e4 | 214 | AliUnicorHN *his = new AliUnicorHN(name,fNdim-1,ax); // result histogram |
215 | AliUnicorHN *eis = new AliUnicorHN(eame,fNdim-1,ax); // temporary storage for errors | |
216 | ||
217 | // sum up the content and errors squared | |
218 | ||
61e4657c | 219 | int k[fgkMaxNdim] = {0}; // old hist multiindex |
220 | int m[fgkMaxNdim] = {0}; // new hist multiindex | |
621688e4 | 221 | for (int i=0; i<GetNbinsX(); i++) { |
222 | OneToMul(i,k); | |
223 | if (k[dim]+1<first) continue; | |
224 | if (k[dim]+1>last) continue; | |
225 | n = 0; | |
226 | for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j]; | |
227 | n = his->MulToOne(m); | |
228 | his->AddBinContent(n+1,GetBinContent(i+1)); | |
229 | eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1)); | |
230 | } | |
231 | ||
232 | // combine content and errors in one histogram | |
233 | ||
234 | for (int i=0; i<his->GetNbinsX(); i++) { | |
235 | his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1))); | |
236 | } | |
237 | ||
238 | his->SetLineColor(this->GetLineColor()); | |
239 | his->SetFillColor(this->GetFillColor()); | |
240 | his->SetMarkerColor(this->GetMarkerColor()); | |
241 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
242 | ||
243 | // some cleanup | |
244 | ||
245 | delete eis; | |
246 | return his; | |
247 | } | |
248 | //============================================================================= | |
28eee19b | 249 | TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first, |
250 | const Int_t * const last) const { | |
251 | ||
621688e4 | 252 | // Project on dimension dim. Use only bins between first[i] and last[i]. |
253 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
254 | // first[i]=0 means from the first bin | |
255 | // last[i]=0 means till the last bin | |
256 | ||
257 | if (dim<0 || dim>fNdim-1) return 0; | |
258 | ||
259 | // calculate the projection; lowest index 0 | |
260 | ||
261 | double *yy = new double[fNbins[dim]]; // value | |
262 | double *ey = new double[fNbins[dim]]; // error | |
263 | for (int i=0; i<fNbins[dim]; i++) yy[i]=0; | |
264 | for (int i=0; i<fNbins[dim]; i++) ey[i]=0; | |
265 | Int_t *k = new Int_t[fNdim]; | |
266 | for (int i=0; i<GetNbinsX(); i++) { | |
267 | OneToMul(i,k); | |
268 | int isgood = 1; // bin within the range? | |
269 | for (int j=0; j<fNdim; j++) { | |
270 | if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;} | |
271 | if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;} | |
272 | } | |
273 | if (isgood) { | |
274 | yy[k[dim]]+=GetBinContent(i+1); | |
275 | ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1); | |
276 | } | |
277 | } | |
278 | ||
279 | // make the projection histogram | |
280 | // use name nam if specified; otherwise generate one | |
281 | ||
282 | TH1D *his; | |
61e4657c | 283 | const char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim); |
621688e4 | 284 | if (fAxis[dim].IsVariableBinSize()) |
285 | his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray()); | |
286 | else | |
287 | his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax()); | |
288 | his->SetXTitle(fAxis[dim].GetTitle()); | |
289 | // or his->GetXaxis()->ImportAttributes(ax); | |
290 | his->Sumw2(); | |
291 | his->SetLineColor(this->GetLineColor()); | |
292 | his->SetFillColor(this->GetFillColor()); | |
293 | his->SetMarkerColor(this->GetMarkerColor()); | |
294 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
295 | for (int i=0; i<his->GetNbinsX(); i++) { | |
296 | his->SetBinContent(i+1,yy[i]); | |
297 | his->SetBinError(i+1,sqrt(ey[i])); | |
298 | } | |
299 | ||
300 | // some cleanup | |
301 | ||
302 | delete [] yy; | |
303 | delete [] ey; | |
304 | delete [] k; | |
305 | // if (name!=nam) delete [] name; | |
306 | ||
307 | return his; | |
308 | } | |
309 | //============================================================================= | |
28eee19b | 310 | TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first, |
311 | const Double_t * const last) const { | |
312 | ||
621688e4 | 313 | // Project on dimension dim. Use only bins between first[i] and last[i]. |
314 | ||
315 | if (dim<0 || dim>fNdim-1) return 0; | |
61e4657c | 316 | Int_t kfirst[fgkMaxNdim] = {0}; |
317 | Int_t klast[fgkMaxNdim] = {0}; | |
621688e4 | 318 | for (int i=0; i<fNdim; i++) { |
28eee19b | 319 | kfirst[i] = fAxis[i].FindFixBin(first[i]); |
320 | klast[i] = fAxis[i].FindFixBin(last[i]); | |
621688e4 | 321 | } |
322 | return ProjectOn(nam,dim,kfirst,klast); | |
323 | } | |
324 | //============================================================================= | |
28eee19b | 325 | TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * |
326 | const first, const Int_t * const last) const { | |
327 | ||
621688e4 | 328 | // Project on dim1 vs dim0. Use only bins between first[i] and last[i]. |
329 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
330 | // first[i]=0 means from the first bin | |
331 | // last[i]=0 means till the last bin | |
332 | ||
333 | if (dim0<0 || dim0>fNdim-1) return 0; | |
334 | if (dim1<0 || dim1>fNdim-1) return 0; | |
335 | ||
336 | // calculate the projection | |
337 | ||
338 | double **yy = new double*[fNbins[dim0]]; // value | |
339 | double **ey = new double*[fNbins[dim0]]; // error | |
340 | for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]]; | |
341 | for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]]; | |
342 | for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0; | |
343 | for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0; | |
344 | Int_t *k = new Int_t[fNdim]; | |
345 | for (int i=0; i<GetNbinsX(); i++) { | |
346 | OneToMul(i,k); | |
347 | int isgood = 1; // bin within the range? | |
348 | for (int j=0; j<fNdim; j++) { | |
349 | if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;} | |
350 | if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;} | |
351 | } | |
352 | if (isgood) { | |
353 | yy[k[dim0]][k[dim1]]+=GetBinContent(i+1); | |
354 | ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1); | |
355 | } | |
356 | } | |
357 | ||
358 | // make the projection histogram | |
359 | // use name nam if specified; otherwise generate one | |
360 | ||
361 | TH2D *his=0; | |
61e4657c | 362 | const char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0); |
621688e4 | 363 | if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) |
364 | his = new TH2D(name,name, | |
365 | fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), | |
366 | fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); | |
367 | else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) | |
368 | his = new TH2D(name,name, | |
369 | fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), | |
370 | fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); | |
371 | else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) | |
372 | his = new TH2D(name,name, | |
373 | fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), | |
374 | fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); | |
375 | else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) | |
376 | his = new TH2D(name,name, | |
377 | fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), | |
378 | fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); | |
379 | his->SetXTitle(fAxis[dim0].GetTitle()); | |
380 | his->SetYTitle(fAxis[dim1].GetTitle()); | |
381 | // or his->GetXaxis()->ImportAttributes(ax0); | |
382 | // or his->GetYaxis()->ImportAttributes(ax1); | |
383 | his->Sumw2(); | |
384 | his->SetLineColor(this->GetLineColor()); | |
385 | his->SetFillColor(this->GetFillColor()); | |
386 | his->SetMarkerColor(this->GetMarkerColor()); | |
387 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
388 | for (int i=0; i<his->GetNbinsX(); i++) | |
389 | for (int j=0; j<his->GetNbinsY(); j++) { | |
390 | his->SetBinContent(i+1,j+1,yy[i][j]); | |
391 | his->SetBinError(i+1,j+1,sqrt(ey[i][j])); | |
392 | } | |
393 | ||
394 | // some cleanup | |
395 | ||
396 | for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i]; | |
397 | for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i]; | |
398 | delete [] yy; | |
399 | delete [] ey; | |
400 | delete [] k; | |
401 | // if (name!=nam) delete [] name; | |
402 | ||
403 | return his; | |
404 | } | |
405 | //============================================================================= | |
28eee19b | 406 | TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Double_t * |
407 | const first, const Double_t * const last) const { | |
408 | ||
409 | // Project on dim1 vs dim0. Use only bins between first[i] and last[i]. | |
410 | ||
411 | if (dim0<0 || dim0>fNdim-1) return 0; | |
412 | if (dim1<0 || dim1>fNdim-1) return 0; | |
413 | Int_t kfirst[fgkMaxNdim] = {0}; | |
414 | Int_t klast[fgkMaxNdim] = {0}; | |
415 | for (int i=0; i<fNdim; i++) { | |
416 | kfirst[i] = fAxis[i].FindFixBin(first[i]); | |
417 | klast[i] = fAxis[i].FindFixBin(last[i]); | |
418 | } | |
419 | return ProjectOn(nam,dim0,dim1,kfirst,klast); | |
420 | } | |
421 | //============================================================================= |