]>
Commit | Line | Data |
---|---|---|
7148817a | 1 | // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007 |
2 | ||
3 | //============================================================================= | |
4 | // multidimensional histogram | |
5 | // Physically, the data are kept in one single one-dimensional histogram. | |
6 | // Projecting on 1, 2, and n-1 dimensions is implemented. | |
7 | // The histogram can be saved on root file as the one-dimensional data | |
8 | // histogram and the axes, thus eternal forward compatibility is ensured. | |
9 | //============================================================================= | |
10 | ||
5a6d201c | 11 | #include <cmath> |
a6e0ebfe | 12 | #include <stdlib.h> |
13 | ||
7148817a | 14 | #include <TFile.h> |
15 | #include <TDirectory.h> | |
16 | #include <TAxis.h> | |
17 | #include <TH2.h> | |
7148817a | 18 | #include "AliDHN.h" |
19 | ||
20 | ClassImp(AliDHN) | |
21 | ||
22 | //============================================================================= | |
23 | AliDHN::AliDHN(Char_t *nam, Int_t ndim, TAxis **ax) | |
24 | : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5), | |
25 | fNdim(ndim) | |
26 | { | |
27 | // Constructor for building from scratch. | |
5a6d201c | 28 | |
7148817a | 29 | // Above, we just have managed to create a single one-dimensional histogram |
5a6d201c | 30 | // with number of bins equal to the product of the numbers of bins in all |
7148817a | 31 | // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5. |
32 | ||
33 | for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]); | |
34 | ||
35 | // for speed reasons, number of bins of each axis is stored on fNbins | |
36 | // and the running product of the latter on fMbins | |
37 | // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1} | |
38 | ||
39 | fMbins[fNdim-1] = 1; | |
40 | for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins(); | |
41 | for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i]; | |
42 | printf("%d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX()); | |
43 | } | |
44 | //============================================================================= | |
45 | AliDHN::AliDHN(Char_t *filnam, Char_t *nam) | |
46 | : TH1D(*((TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"))), | |
47 | fNdim(0) | |
48 | { | |
49 | // Constructor for reading from file. | |
50 | ||
51 | TFile *f = TFile::Open(filnam,"read"); | |
52 | f->cd(nam); | |
53 | TAxis *ax[fMaxNdim]; | |
54 | char axnam[1000]; | |
55 | for (fNdim=0; fNdim<fMaxNdim; fNdim++) { | |
56 | sprintf(axnam,"axis%d",fNdim); | |
57 | ax[fNdim] = (TAxis *) gDirectory->Get(axnam); | |
58 | if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]); | |
59 | else break; | |
60 | } | |
61 | f->Close(); | |
62 | ||
63 | fMbins[fNdim-1] = 1; | |
64 | for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins(); | |
65 | for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i]; | |
66 | ||
67 | if (GetNbinsX()!=Albins(fNdim,ax)) { | |
68 | printf("number of bins of histo %d differs from product of nbins of axes %d\n", | |
69 | GetNbinsX(),Albins(fNdim,ax)); | |
70 | printf("bombing\n"); | |
71 | exit(-1); | |
72 | } | |
73 | printf("%d-dimensional histogram %s with %d bins read from file %s\n", | |
74 | fNdim,nam,GetNbinsX(),filnam); | |
75 | } | |
76 | //============================================================================= | |
77 | Int_t AliDHN::Albins(Int_t n, TAxis **ax) | |
78 | { | |
79 | // Calculate product of nbins of ax[0]...ax[n-1] | |
80 | // (= total number of bins of the multidimensional histogram to be made). | |
81 | ||
82 | Int_t nbins = 1; | |
83 | // while (n--) nbins *= ax[n]->GetNbins(); | |
84 | for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins(); | |
85 | return nbins; | |
86 | } | |
87 | //============================================================================= | |
88 | Int_t AliDHN::MulToOne(Int_t *k) const | |
89 | { | |
90 | // Calculate the 1-dim index n from n-dim indices k[fNdim]. | |
91 | // Valid k[i] should be between 0 and fNbins[i]-1. | |
92 | // Valid result will be between 0 and GetNbinsX()-1. | |
93 | // Return -1 if under- or overflow in any dimension. | |
94 | ||
95 | Int_t n = 0; | |
96 | for (int i=0; i<fNdim; i++) { | |
97 | if (k[i]<0) return -1; | |
98 | if (k[i]>=fNbins[i]) return -1; | |
99 | n += fMbins[i]*k[i]; | |
100 | } | |
101 | return n; | |
102 | } | |
103 | //============================================================================= | |
104 | Int_t AliDHN::MulToOne(Double_t *x) | |
105 | { | |
106 | // Calculate the 1-dim index n from n-dim vector x, representing the | |
107 | // abscissa of the n-dim histogram. The result will be between 0 and | |
108 | // GetNbinsX()-1. | |
109 | ||
110 | Int_t k[fMaxNdim]; | |
111 | for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1; | |
112 | return MulToOne(k); | |
113 | } | |
114 | //============================================================================= | |
5a6d201c | 115 | void AliDHN::OneToMul(Int_t n, Int_t *k) const |
7148817a | 116 | { |
117 | // Calculate the n-dim indices k[fNdim] from 1-dim index n. | |
118 | // Valid n should be between 0 and GetNbinsX()-1. | |
119 | // Valid results will be between 0 and fNbins[i]-1. | |
120 | ||
121 | div_t idi; // integer division structure | |
122 | for (int i=0; i<fNdim; i++) { | |
123 | idi = div(n,fMbins[i]); | |
124 | k[i] = idi.quot; // quotient | |
125 | n = idi.rem; // remainder | |
126 | } | |
127 | } | |
128 | //============================================================================= | |
5a6d201c | 129 | Int_t AliDHN::Fill(Double_t *xx, Double_t w) |
7148817a | 130 | { |
131 | // Fill the histogram. The array xx holds the abscissa information, w is the | |
132 | // weigth. The 1-dim histogram is filled using the standard Fill method | |
133 | // (direct access to the arrays was tried and was not faster). | |
134 | ||
135 | int nbin = MulToOne(xx); | |
5a6d201c | 136 | if (nbin == -1) return 0; |
137 | return TH1D::Fill(nbin+1,w); | |
7148817a | 138 | } |
139 | //============================================================================= | |
5a6d201c | 140 | Int_t AliDHN::Fill(Double_t x0, Double_t x1, ...) |
141 | { | |
142 | // Fill the histogram. Arguments are passed as doubles rather than array. | |
143 | va_list ap; | |
144 | Double_t xx[fMaxNdim] = {x0, x1}; | |
145 | va_start(ap,x1); | |
146 | for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t); | |
147 | Double_t weigth = va_arg(ap,Double_t); | |
148 | va_end(ap); | |
149 | return Fill(xx,weigth); | |
150 | } | |
151 | //============================================================================= | |
152 | Int_t AliDHN::Write() const | |
7148817a | 153 | { |
154 | // Save the 1-dim histo and the axes in a subdirectory on file. This might | |
155 | // not be the most elegant way but it is very simple and backward and forward | |
156 | // compatible. | |
157 | ||
158 | Int_t nbytes = 0; | |
159 | TH1D histo(*this); | |
160 | histo.SetName("histo"); // hadd bug fix | |
161 | TDirectory *dest = gDirectory->mkdir(GetName()); | |
162 | if (dest) { | |
163 | dest->cd(); | |
164 | nbytes += histo.Write("histo"); | |
165 | char axnam[1000]; | |
166 | for (int i=0; i<fNdim; i++) { | |
167 | sprintf(axnam,"axis%d",i); | |
168 | nbytes += fAxis[i].Write(axnam); | |
169 | } | |
170 | printf(" %s stored in %s\n",GetName(),dest->GetPath()); | |
171 | dest->cd(".."); | |
172 | } | |
173 | return nbytes; | |
174 | } | |
175 | //============================================================================= | |
176 | AliDHN *AliDHN::ProjectAlong(char *nam, Int_t dim, Int_t first, Int_t last) | |
177 | { | |
178 | // Reduce dimension dim by summing up its bins between first and last. | |
179 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
180 | // Return the resulting fNdim-1 dimensional histogram. | |
181 | ||
182 | if (dim<0 || dim>fNdim-1) return 0; | |
183 | if (first<0) first = 1; | |
184 | if (last<0) last = fNbins[dim]; | |
185 | ||
186 | // create new (reduced) histogram | |
187 | ||
188 | TAxis *ax[fMaxNdim-1]; | |
189 | int n=0; | |
190 | for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i); | |
191 | char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim); | |
192 | char *eame = Form("%s_error",name); | |
193 | AliDHN *his = new AliDHN(name,fNdim-1,ax); // result histogram | |
194 | AliDHN *eis = new AliDHN(eame,fNdim-1,ax); // temporary storage for errors | |
195 | ||
196 | // sum up the content and errors squared | |
197 | ||
198 | Int_t *k = new Int_t[fNdim]; // old hist multiindex | |
199 | Int_t *m = new Int_t[fNdim-1]; // new hist multiindex | |
200 | for (int i=0; i<GetNbinsX(); i++) { | |
201 | OneToMul(i,k); | |
202 | if (k[dim]+1<first) continue; | |
203 | if (k[dim]+1>last) continue; | |
5a6d201c | 204 | n = 0; |
7148817a | 205 | for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j]; |
206 | n = his->MulToOne(m); | |
207 | his->AddBinContent(n+1,GetBinContent(i+1)); | |
208 | eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1)); | |
209 | } | |
210 | ||
211 | // combine content and errors in one histogram | |
212 | ||
213 | for (int i=0; i<his->GetNbinsX(); i++) { | |
5a6d201c | 214 | his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1))); |
7148817a | 215 | } |
216 | ||
217 | his->SetLineColor(this->GetLineColor()); | |
218 | his->SetFillColor(this->GetFillColor()); | |
219 | his->SetMarkerColor(this->GetMarkerColor()); | |
220 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
221 | ||
222 | // some cleanup | |
223 | ||
224 | delete eis; | |
225 | return his; | |
226 | } | |
227 | //============================================================================= | |
5a6d201c | 228 | TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Int_t *first, Int_t *last) const |
7148817a | 229 | { |
230 | // Project on dimension dim. Use only bins between first[i] and last[i]. | |
231 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
232 | // first[i]=0 means from the first bin | |
233 | // last[i]=0 means till the last bin | |
234 | ||
235 | if (dim<0 || dim>fNdim-1) return 0; | |
236 | ||
237 | // calculate the projection; lowest index 0 | |
238 | ||
239 | double *yy = new double[fNbins[dim]]; // value | |
240 | double *ey = new double[fNbins[dim]]; // error | |
241 | for (int i=0; i<fNbins[dim]; i++) yy[i]=0; | |
242 | for (int i=0; i<fNbins[dim]; i++) ey[i]=0; | |
243 | Int_t *k = new Int_t[fNdim]; | |
244 | for (int i=0; i<GetNbinsX(); i++) { | |
245 | OneToMul(i,k); | |
246 | int isgood = 1; // bin within the range? | |
247 | for (int j=0; j<fNdim; j++) { | |
248 | if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;} | |
249 | if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;} | |
250 | } | |
251 | if (isgood) { | |
252 | yy[k[dim]]+=GetBinContent(i+1); | |
253 | ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1); | |
254 | } | |
255 | } | |
256 | ||
257 | // make the projection histogram | |
258 | // use name nam if specified; otherwise generate one | |
259 | ||
260 | TH1D *his; | |
261 | char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim); | |
262 | if (fAxis[dim].IsVariableBinSize()) | |
263 | his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray()); | |
264 | else | |
265 | his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax()); | |
266 | his->SetXTitle(fAxis[dim].GetTitle()); | |
267 | // or his->GetXaxis()->ImportAttributes(ax); | |
268 | his->Sumw2(); | |
269 | his->SetLineColor(this->GetLineColor()); | |
270 | his->SetFillColor(this->GetFillColor()); | |
271 | his->SetMarkerColor(this->GetMarkerColor()); | |
272 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
273 | for (int i=0; i<his->GetNbinsX(); i++) { | |
274 | his->SetBinContent(i+1,yy[i]); | |
5a6d201c | 275 | his->SetBinError(i+1,sqrt(ey[i])); |
7148817a | 276 | } |
277 | ||
278 | // some cleanup | |
279 | ||
280 | delete [] yy; | |
281 | delete [] ey; | |
282 | delete [] k; | |
283 | // if (name!=nam) delete [] name; | |
284 | ||
285 | return his; | |
286 | } | |
287 | //============================================================================= | |
288 | TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Double_t *first, Double_t *last) | |
289 | { | |
290 | // Project on dimension dim. Use only bins between first[i] and last[i]. | |
291 | ||
292 | if (dim<0 || dim>fNdim-1) return 0; | |
293 | Int_t kfirst[fMaxNdim]; | |
294 | Int_t klast[fMaxNdim]; | |
295 | for (int i=0; i<fNdim; i++) { | |
296 | kfirst[i] = fAxis[i].FindBin(first[i]); | |
297 | klast[i] = fAxis[i].FindBin(last[i]); | |
298 | } | |
299 | return ProjectOn(nam,dim,kfirst,klast); | |
300 | } | |
301 | //============================================================================= | |
5a6d201c | 302 | TH2D *AliDHN::ProjectOn(char *nam, Int_t dim0, Int_t dim1, Int_t *first, Int_t *last) const |
7148817a | 303 | { |
304 | // Project on dim1 vs dim0. Use only bins between first[i] and last[i]. | |
305 | // Use root convention: bin=1 is the first bin, bin=nbins is the last. | |
306 | // first[i]=0 means from the first bin | |
307 | // last[i]=0 means till the last bin | |
308 | ||
309 | if (dim0<0 || dim0>fNdim-1) return 0; | |
310 | if (dim1<0 || dim1>fNdim-1) return 0; | |
311 | ||
312 | // calculate the projection | |
313 | ||
314 | double **yy = new double*[fNbins[dim0]]; // value | |
315 | double **ey = new double*[fNbins[dim0]]; // error | |
316 | for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]]; | |
317 | for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]]; | |
318 | for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0; | |
319 | for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0; | |
320 | Int_t *k = new Int_t[fNdim]; | |
321 | for (int i=0; i<GetNbinsX(); i++) { | |
322 | OneToMul(i,k); | |
323 | int isgood = 1; // bin within the range? | |
324 | for (int j=0; j<fNdim; j++) { | |
325 | if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;} | |
326 | if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;} | |
327 | } | |
328 | if (isgood) { | |
329 | yy[k[dim0]][k[dim1]]+=GetBinContent(i+1); | |
330 | ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1); | |
331 | } | |
332 | } | |
333 | ||
334 | // make the projection histogram | |
335 | // use name nam if specified; otherwise generate one | |
336 | ||
337 | TH2D *his=0; | |
338 | char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0); | |
339 | if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) | |
340 | his = new TH2D(name,name, | |
341 | fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), | |
342 | fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); | |
343 | else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) | |
344 | his = new TH2D(name,name, | |
345 | fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), | |
346 | fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); | |
347 | else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) | |
348 | his = new TH2D(name,name, | |
349 | fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), | |
350 | fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); | |
351 | else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) | |
352 | his = new TH2D(name,name, | |
353 | fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), | |
354 | fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); | |
355 | his->SetXTitle(fAxis[dim0].GetTitle()); | |
356 | his->SetYTitle(fAxis[dim1].GetTitle()); | |
357 | // or his->GetXaxis()->ImportAttributes(ax0); | |
358 | // or his->GetYaxis()->ImportAttributes(ax1); | |
359 | his->Sumw2(); | |
360 | his->SetLineColor(this->GetLineColor()); | |
361 | his->SetFillColor(this->GetFillColor()); | |
362 | his->SetMarkerColor(this->GetMarkerColor()); | |
363 | his->SetMarkerStyle(this->GetMarkerStyle()); | |
364 | for (int i=0; i<his->GetNbinsX(); i++) | |
365 | for (int j=0; j<his->GetNbinsY(); j++) { | |
366 | his->SetBinContent(i+1,j+1,yy[i][j]); | |
5a6d201c | 367 | his->SetBinError(i+1,j+1,sqrt(ey[i][j])); |
7148817a | 368 | } |
369 | ||
370 | // some cleanup | |
371 | ||
372 | for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i]; | |
373 | for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i]; | |
374 | delete [] yy; | |
375 | delete [] ey; | |
376 | delete [] k; | |
377 | // if (name!=nam) delete [] name; | |
378 | ||
379 | return his; | |
380 | } | |
381 | //============================================================================= |