]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************** | |
2 | * Copyright(c) 1998-1999, 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 | /* $Id: AliTHn.cxx 20164 2007-08-14 15:31:50Z morsch $ */ | |
17 | ||
18 | // Use AliTHn instead of AliCFContainer and your memory consumption will be drastically reduced | |
19 | // As AliTHn derives from AliCFContainer, you can just replace your current AliCFContainer object by AliTHn | |
20 | // Once you have the merged output, call FillParent() and you can use AliCFContainer as usual | |
21 | // | |
22 | // this storage container is optimized for small memory usage | |
23 | // under/over flow bins do not exist | |
24 | // sumw2 structure is float only and only create when the weight != 1 | |
25 | // all histogram functionality (projections, axis, ranges, etc) are taken from THnSparse by propagating | |
26 | // the information up into the THnSparse structure (in the ananalysis *after* data processing and merging) | |
27 | // | |
28 | // the derivation from THnSparse is obviously against many OO rules. correct would be a common baseclass of THnSparse and THn. | |
29 | // | |
30 | // Author: Jan Fiete Grosse-Oetringhaus | |
31 | ||
32 | #include "AliTHn.h" | |
33 | #include "TList.h" | |
34 | #include "TCollection.h" | |
35 | #include "AliLog.h" | |
36 | #include "TArrayF.h" | |
37 | #include "THnSparse.h" | |
38 | #include "TMath.h" | |
39 | ||
40 | ClassImp(AliTHn) | |
41 | ||
42 | AliTHn::AliTHn() : | |
43 | AliCFContainer(), | |
44 | fNBins(0), | |
45 | fNVars(0), | |
46 | fNSteps(0), | |
47 | fValues(0), | |
48 | fSumw2(0), | |
49 | axisCache(0), | |
50 | fNbinsCache(0), | |
51 | fLastVars(0), | |
52 | fLastBins(0) | |
53 | { | |
54 | // Constructor | |
55 | } | |
56 | ||
57 | AliTHn::AliTHn(const Char_t* name, const Char_t* title,const Int_t nSelStep, const Int_t nVarIn, const Int_t* nBinIn) : | |
58 | AliCFContainer(name, title, nSelStep, nVarIn, nBinIn), | |
59 | fNBins(0), | |
60 | fNVars(nVarIn), | |
61 | fNSteps(nSelStep), | |
62 | fValues(0), | |
63 | fSumw2(0), | |
64 | axisCache(0), | |
65 | fNbinsCache(0), | |
66 | fLastVars(0), | |
67 | fLastBins(0) | |
68 | { | |
69 | // Constructor | |
70 | ||
71 | fNBins = 1; | |
72 | for (Int_t i=0; i<fNVars; i++) | |
73 | fNBins *= nBinIn[i]; | |
74 | ||
75 | Init(); | |
76 | } | |
77 | ||
78 | void AliTHn::Init() | |
79 | { | |
80 | // initialize | |
81 | ||
82 | fValues = new TArrayF*[fNSteps]; | |
83 | fSumw2 = new TArrayF*[fNSteps]; | |
84 | ||
85 | for (Int_t i=0; i<fNSteps; i++) | |
86 | { | |
87 | fValues[i] = 0; | |
88 | fSumw2[i] = 0; | |
89 | } | |
90 | } | |
91 | ||
92 | AliTHn::AliTHn(const AliTHn &c) : | |
93 | AliCFContainer(c), | |
94 | fNBins(c.fNBins), | |
95 | fNVars(c.fNVars), | |
96 | fNSteps(c.fNSteps), | |
97 | fValues(new TArrayF*[c.fNSteps]), | |
98 | fSumw2(new TArrayF*[c.fNSteps]), | |
99 | axisCache(0), | |
100 | fNbinsCache(0), | |
101 | fLastVars(0), | |
102 | fLastBins(0) | |
103 | { | |
104 | // | |
105 | // AliTHn copy constructor | |
106 | // | |
107 | ||
108 | memset(fValues,0,fNSteps*sizeof(TArrayF*)); | |
109 | memset(fSumw2,0,fNSteps*sizeof(TArrayF*)); | |
110 | ||
111 | for (Int_t i=0; i<fNSteps; i++) { | |
112 | if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i])); | |
113 | if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i])); | |
114 | } | |
115 | ||
116 | } | |
117 | ||
118 | AliTHn::~AliTHn() | |
119 | { | |
120 | // Destructor | |
121 | ||
122 | DeleteContainers(); | |
123 | ||
124 | delete[] fValues; | |
125 | delete[] fSumw2; | |
126 | delete[] axisCache; | |
127 | delete[] fNbinsCache; | |
128 | delete[] fLastVars; | |
129 | delete[] fLastBins; | |
130 | } | |
131 | ||
132 | void AliTHn::DeleteContainers() | |
133 | { | |
134 | // delete data containers | |
135 | ||
136 | for (Int_t i=0; i<fNSteps; i++) | |
137 | { | |
138 | if (fValues && fValues[i]) | |
139 | { | |
140 | delete fValues[i]; | |
141 | fValues[i] = 0; | |
142 | } | |
143 | ||
144 | if (fSumw2 && fSumw2[i]) | |
145 | { | |
146 | delete fSumw2[i]; | |
147 | fSumw2[i] = 0; | |
148 | } | |
149 | } | |
150 | } | |
151 | ||
152 | //____________________________________________________________________ | |
153 | AliTHn &AliTHn::operator=(const AliTHn &c) | |
154 | { | |
155 | // assigment operator | |
156 | ||
157 | if (this != &c) { | |
158 | AliCFContainer::operator=(c); | |
159 | fNBins=c.fNBins; | |
160 | fNVars=c.fNVars; | |
161 | if(fNSteps) { | |
162 | for(Int_t i=0; i< fNSteps; ++i) { | |
163 | delete fValues[i]; | |
164 | delete fSumw2[i]; | |
165 | } | |
166 | delete [] fValues; | |
167 | delete [] fSumw2; | |
168 | } | |
169 | fNSteps=c.fNSteps; | |
170 | if(fNSteps) { | |
171 | fValues=new TArrayF*[fNSteps]; | |
172 | fSumw2=new TArrayF*[fNSteps]; | |
173 | memset(fValues,0,fNSteps*sizeof(TArrayF*)); | |
174 | memset(fSumw2,0,fNSteps*sizeof(TArrayF*)); | |
175 | ||
176 | for (Int_t i=0; i<fNSteps; i++) { | |
177 | if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i])); | |
178 | if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i])); | |
179 | } | |
180 | } else { | |
181 | fValues = 0; | |
182 | fSumw2 = 0; | |
183 | } | |
184 | delete [] axisCache; | |
185 | axisCache = new TAxis*[fNVars]; | |
186 | memcpy(axisCache, c.axisCache, fNVars*sizeof(TAxis*)); | |
187 | } | |
188 | return *this; | |
189 | } | |
190 | ||
191 | //____________________________________________________________________ | |
192 | void AliTHn::Copy(TObject& c) const | |
193 | { | |
194 | // copy function | |
195 | ||
196 | AliTHn& target = (AliTHn &) c; | |
197 | ||
198 | AliCFContainer::Copy(target); | |
199 | ||
200 | target.fNSteps = fNSteps; | |
201 | target.fNBins = fNBins; | |
202 | target.fNVars = fNVars; | |
203 | ||
204 | target.Init(); | |
205 | ||
206 | for (Int_t i=0; i<fNSteps; i++) | |
207 | { | |
208 | if (fValues[i]) | |
209 | target.fValues[i] = new TArrayF(*(fValues[i])); | |
210 | else | |
211 | target.fValues[i] = 0; | |
212 | ||
213 | if (fSumw2[i]) | |
214 | target.fSumw2[i] = new TArrayF(*(fSumw2[i])); | |
215 | else | |
216 | target.fSumw2[i] = 0; | |
217 | } | |
218 | } | |
219 | ||
220 | //____________________________________________________________________ | |
221 | Long64_t AliTHn::Merge(TCollection* list) | |
222 | { | |
223 | // Merge a list of AliTHn objects with this (needed for | |
224 | // PROOF). | |
225 | // Returns the number of merged objects (including this). | |
226 | ||
227 | if (!list) | |
228 | return 0; | |
229 | ||
230 | if (list->IsEmpty()) | |
231 | return 1; | |
232 | ||
233 | AliCFContainer::Merge(list); | |
234 | ||
235 | TIterator* iter = list->MakeIterator(); | |
236 | TObject* obj; | |
237 | ||
238 | Int_t count = 0; | |
239 | while ((obj = iter->Next())) { | |
240 | ||
241 | AliTHn* entry = dynamic_cast<AliTHn*> (obj); | |
242 | if (entry == 0) | |
243 | continue; | |
244 | ||
245 | for (Int_t i=0; i<fNSteps; i++) | |
246 | { | |
247 | if (entry->fValues[i]) | |
248 | { | |
249 | if (!fValues[i]) | |
250 | fValues[i] = new TArrayF(fNBins); | |
251 | ||
252 | for (Long64_t l = 0; l<fNBins; l++) | |
253 | fValues[i]->GetArray()[l] += entry->fValues[i]->GetArray()[l]; | |
254 | } | |
255 | ||
256 | if (entry->fSumw2[i]) | |
257 | { | |
258 | if (!fSumw2[i]) | |
259 | fSumw2[i] = new TArrayF(fNBins); | |
260 | ||
261 | for (Long64_t l = 0; l<fNBins; l++) | |
262 | fSumw2[i]->GetArray()[l] += entry->fSumw2[i]->GetArray()[l]; | |
263 | } | |
264 | } | |
265 | ||
266 | count++; | |
267 | } | |
268 | ||
269 | return count+1; | |
270 | } | |
271 | ||
272 | void AliTHn::Fill(const Double_t *var, Int_t istep, Double_t weight) | |
273 | { | |
274 | // fills an entry | |
275 | ||
276 | // fill axis cache | |
277 | if (!axisCache) | |
278 | { | |
279 | axisCache = new TAxis*[fNVars]; | |
280 | fNbinsCache = new Int_t[fNVars]; | |
281 | for (Int_t i=0; i<fNVars; i++) | |
282 | { | |
283 | axisCache[i] = GetAxis(i, 0); | |
284 | fNbinsCache[i] = axisCache[i]->GetNbins(); | |
285 | } | |
286 | ||
287 | fLastVars = new Double_t[fNVars]; | |
288 | fLastBins = new Int_t[fNVars]; | |
289 | ||
290 | // initial values to prevent checking for 0 below | |
291 | for (Int_t i=0; i<fNVars; i++) | |
292 | { | |
293 | fLastBins[i] = axisCache[i]->FindBin(var[i]); | |
294 | fLastVars[i] = var[i]; | |
295 | } | |
296 | } | |
297 | ||
298 | // calculate global bin index | |
299 | Long64_t bin = 0; | |
300 | for (Int_t i=0; i<fNVars; i++) | |
301 | { | |
302 | bin *= fNbinsCache[i]; | |
303 | ||
304 | Int_t tmpBin = 0; | |
305 | if (fLastVars[i] == var[i]) | |
306 | tmpBin = fLastBins[i]; | |
307 | else | |
308 | { | |
309 | tmpBin = axisCache[i]->FindBin(var[i]); | |
310 | fLastBins[i] = tmpBin; | |
311 | fLastVars[i] = var[i]; | |
312 | } | |
313 | //Printf("%d", tmpBin); | |
314 | ||
315 | // under/overflow not supported | |
316 | if (tmpBin < 1 || tmpBin > fNbinsCache[i]) | |
317 | return; | |
318 | ||
319 | // bins start from 0 here | |
320 | bin += tmpBin - 1; | |
321 | // Printf("%lld", bin); | |
322 | } | |
323 | ||
324 | if (!fValues[istep]) | |
325 | { | |
326 | fValues[istep] = new TArrayF(fNBins); | |
327 | AliInfo(Form("Created values container for step %d", istep)); | |
328 | } | |
329 | ||
330 | if (weight != 1) | |
331 | { | |
332 | // initialize with already filled entries (which have been filled with weight == 1), in this case fSumw2 := fValues | |
333 | if (!fSumw2[istep]) | |
334 | { | |
335 | fSumw2[istep] = new TArrayF(*fValues[istep]); | |
336 | AliInfo(Form("Created sumw2 container for step %d", istep)); | |
337 | } | |
338 | } | |
339 | ||
340 | fValues[istep]->GetArray()[bin] += weight; | |
341 | if (fSumw2[istep]) | |
342 | fSumw2[istep]->GetArray()[bin] += weight * weight; | |
343 | ||
344 | // Printf("%f", fValues[istep][bin]); | |
345 | ||
346 | // debug | |
347 | // AliCFContainer::Fill(var, istep, weight); | |
348 | } | |
349 | ||
350 | Long64_t AliTHn::GetGlobalBinIndex(const Int_t* binIdx) | |
351 | { | |
352 | // calculates global bin index | |
353 | // binIdx contains TAxis bin indexes | |
354 | // here bin count starts at 0 because we do not have over/underflow bins | |
355 | ||
356 | Long64_t bin = 0; | |
357 | for (Int_t i=0; i<fNVars; i++) | |
358 | { | |
359 | bin *= GetAxis(i, 0)->GetNbins(); | |
360 | bin += binIdx[i] - 1; | |
361 | } | |
362 | ||
363 | return bin; | |
364 | } | |
365 | ||
366 | void AliTHn::FillContainer(AliCFContainer* cont) | |
367 | { | |
368 | // fills the information stored in the buffer in this class into the container <cont> | |
369 | ||
370 | for (Int_t i=0; i<fNSteps; i++) | |
371 | { | |
372 | if (!fValues[i]) | |
373 | continue; | |
374 | ||
375 | Float_t* source = fValues[i]->GetArray(); | |
376 | // if fSumw2 is not stored, the sqrt of the number of bin entries in source is filled below; otherwise we use fSumw2 | |
377 | Float_t* sourceSumw2 = source; | |
378 | if (fSumw2[i]) | |
379 | sourceSumw2 = fSumw2[i]->GetArray(); | |
380 | ||
381 | THnSparse* target = cont->GetGrid(i)->GetGrid(); | |
382 | ||
383 | Int_t* binIdx = new Int_t[fNVars]; | |
384 | Int_t* nBins = new Int_t[fNVars]; | |
385 | for (Int_t j=0; j<fNVars; j++) | |
386 | { | |
387 | binIdx[j] = 1; | |
388 | nBins[j] = target->GetAxis(j)->GetNbins(); | |
389 | } | |
390 | ||
391 | Long64_t count = 0; | |
392 | ||
393 | while (1) | |
394 | { | |
395 | // for (Int_t j=0; j<fNVars; j++) | |
396 | // printf("%d ", binIdx[j]); | |
397 | ||
398 | Long64_t globalBin = GetGlobalBinIndex(binIdx); | |
399 | // Printf(" --> %lld", globalBin); | |
400 | ||
401 | if (source[globalBin] != 0) | |
402 | { | |
403 | target->SetBinContent(binIdx, source[globalBin]); | |
404 | target->SetBinError(binIdx, TMath::Sqrt(sourceSumw2[globalBin])); | |
405 | ||
406 | count++; | |
407 | } | |
408 | ||
409 | binIdx[fNVars-1]++; | |
410 | ||
411 | for (Int_t j=fNVars-1; j>0; j--) | |
412 | { | |
413 | if (binIdx[j] > nBins[j]) | |
414 | { | |
415 | binIdx[j] = 1; | |
416 | binIdx[j-1]++; | |
417 | } | |
418 | } | |
419 | ||
420 | if (binIdx[0] > nBins[0]) | |
421 | break; | |
422 | } | |
423 | ||
424 | AliInfo(Form("Step %d: copied %lld entries out of %lld bins", i, count, GetGlobalBinIndex(binIdx))); | |
425 | ||
426 | delete[] binIdx; | |
427 | delete[] nBins; | |
428 | } | |
429 | } | |
430 | ||
431 | void AliTHn::FillParent() | |
432 | { | |
433 | // fills the information stored in the buffer in this class into the baseclass containers | |
434 | ||
435 | FillContainer(this); | |
436 | } | |
437 | ||
438 | void AliTHn::ReduceAxis() | |
439 | { | |
440 | // "removes" one axis by summing over the axis and putting the entry to bin 1 | |
441 | // TODO presently only implemented for the last axis | |
442 | ||
443 | Int_t axis = fNVars-1; | |
444 | ||
445 | for (Int_t i=0; i<fNSteps; i++) | |
446 | { | |
447 | if (!fValues[i]) | |
448 | continue; | |
449 | ||
450 | Float_t* source = fValues[i]->GetArray(); | |
451 | Float_t* sourceSumw2 = 0; | |
452 | if (fSumw2[i]) | |
453 | sourceSumw2 = fSumw2[i]->GetArray(); | |
454 | ||
455 | THnSparse* target = GetGrid(i)->GetGrid(); | |
456 | ||
457 | Int_t* binIdx = new Int_t[fNVars]; | |
458 | Int_t* nBins = new Int_t[fNVars]; | |
459 | for (Int_t j=0; j<fNVars; j++) | |
460 | { | |
461 | binIdx[j] = 1; | |
462 | nBins[j] = target->GetAxis(j)->GetNbins(); | |
463 | } | |
464 | ||
465 | Long64_t count = 0; | |
466 | ||
467 | while (1) | |
468 | { | |
469 | // sum over axis <axis> | |
470 | Float_t sumValues = 0; | |
471 | Float_t sumSumw2 = 0; | |
472 | for (Int_t j=1; j<=nBins[axis]; j++) | |
473 | { | |
474 | binIdx[axis] = j; | |
475 | Long64_t globalBin = GetGlobalBinIndex(binIdx); | |
476 | sumValues += source[globalBin]; | |
477 | source[globalBin] = 0; | |
478 | ||
479 | if (sourceSumw2) | |
480 | { | |
481 | sumSumw2 += sourceSumw2[globalBin]; | |
482 | sourceSumw2[globalBin] = 0; | |
483 | } | |
484 | } | |
485 | binIdx[axis] = 1; | |
486 | ||
487 | Long64_t globalBin = GetGlobalBinIndex(binIdx); | |
488 | source[globalBin] = sumValues; | |
489 | if (sourceSumw2) | |
490 | sourceSumw2[globalBin] = sumSumw2; | |
491 | ||
492 | count++; | |
493 | ||
494 | // next bin | |
495 | binIdx[fNVars-2]++; | |
496 | ||
497 | for (Int_t j=fNVars-2; j>0; j--) | |
498 | { | |
499 | if (binIdx[j] > nBins[j]) | |
500 | { | |
501 | binIdx[j] = 1; | |
502 | binIdx[j-1]++; | |
503 | } | |
504 | } | |
505 | ||
506 | if (binIdx[0] > nBins[0]) | |
507 | break; | |
508 | } | |
509 | ||
510 | AliInfo(Form("Step %d: reduced %lld bins to %lld entries", i, GetGlobalBinIndex(binIdx), count)); | |
511 | ||
512 | delete[] binIdx; | |
513 | delete[] nBins; | |
514 | } | |
515 | } |