PWG/Tools/AliTHn.cxx

   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 }