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