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