CORRFW/AliCFGridSparse.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 //                                                                    //
  17 // AliCFGridSparse Class                                              //
  18 // Class to accumulate data on an N-dimensional grid, to be used      //
  19 // as input to get corrections for Reconstruction & Trigger efficiency//
  20 // Based on root THnSparse                                            //
  21 // -- Author : S.Arcelli                                              //
  22 // Still to be done:                                                  //
  23 // --Interpolate among bins in a range                                //
  24 //--------------------------------------------------------------------//
  25 //
  26 //
  27 #include "AliCFGridSparse.h"
  28 #include "THnSparse.h"
  29 #include "AliLog.h"
  30 #include "TMath.h"
  31 #include "TROOT.h"
  32 #include "TH1D.h"
  33 #include "TH2D.h"
  34 #include "TH3D.h"
  35 #include "TAxis.h"
  36
  37 //____________________________________________________________________
  38 ClassImp(AliCFGridSparse)
  39
  40 //____________________________________________________________________
  41 AliCFGridSparse::AliCFGridSparse() :
  42   AliCFVGrid(),
  43   fData(0x0)
  44 {
  45   // default constructor
  46 }
  47 //____________________________________________________________________
  48 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
  49   AliCFVGrid(name,title),
  50   fData(0x0)
  51 {
  52   // default constructor
  53 }
  54 //____________________________________________________________________
  55 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, const Int_t nVarIn, const Int_t * nBinIn, const Double_t *binLimitsIn) :
  56   AliCFVGrid(name,title,nVarIn,nBinIn,binLimitsIn),
  57   fData(0x0)
  58 {
  59   //
  60   // main constructor
  61   //
  62
  63   fData=new THnSparseF(name,title,fNVar,fNVarBins);
  64
  65   if(binLimitsIn){
  66     for(Int_t ivar=0;ivar<fNVar;ivar++){
  67       Int_t nbins=fNVarBins[ivar]+1;
  68       Double_t *array= new Double_t[nbins];
  69       for(Int_t i=0;i<nbins;i++){
  70         array[i]=fVarBinLimits[fOffset[ivar]+i];
  71       }
  72       fData->SetBinEdges(ivar, array);
  73       delete array;
  74     }
  75   }
  76 }
  77 //____________________________________________________________________
  78 AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
  79   AliCFVGrid(c),
  80   fData(c.fData)
  81 {
  82   //
  83   // copy constructor
  84   //
  85   ((AliCFGridSparse &)c).Copy(*this);
  86 }
  87
  88 //____________________________________________________________________
  89 AliCFGridSparse::~AliCFGridSparse()
  90 {
  91   //
  92   // destructor
  93   //
  94   if(fData) delete fData;
  95 }
  96
  97 //____________________________________________________________________
  98 AliCFGridSparse &AliCFGridSparse::operator=(const AliCFGridSparse &c)
  99 {
 100   //
 101   // assigment operator
 102   //
 103   if (this != &c)
 104     ((AliCFGridSparse &) c).Copy(*this);
 105   return *this;
 106 }
 107
 108 //____________________________________________________________________
 109 void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t *array)
 110 {
 111   //
 112   // setting the arrays containing the bin limits
 113   //
 114   fData->SetBinEdges(ivar, array);
 115   //then fill the appropriate array in ALICFFrame, to be able to use
 116   //the getter, in case....
 117   Int_t nbins=fNVarBins[ivar]+1;
 118   for(Int_t i=0;i<nbins;i++){
 119     fVarBinLimits[fOffset[ivar]+i] =array[i];
 120   }
 121 }
 122
 123 //____________________________________________________________________
 124 void AliCFGridSparse::Fill(Double_t *var, Double_t weight)
 125 {
 126   //
 127   // Fill the grid,
 128   // given a set of values of the input variable,
 129   // with weight (by default w=1)
 130   //
 131   fData->Fill(var,weight);
 132 }
 133
 134 //___________________________________________________________________
 135 TH1D *AliCFGridSparse::Project(Int_t ivar) const
 136 {
 137   //
 138   // Make a 1D projection along variable ivar
 139   //
 140
 141   //exclude overflows in hidden dimesniosn by default (used in projections)
 142   if(fExclOffEntriesInProj){
 143     for(Int_t i=0;i<fNVar;i++){
 144       if(i!=ivar)
 145         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
 146     }
 147   }
 148   TH1D *hist=fData->Projection(ivar);
 149   return hist;
 150 }
 151 //___________________________________________________________________
 152 TH2D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2) const
 153 {
 154   //
 155   // Make a 2D projection along variables ivar1 & ivar2
 156   //
 157
 158   //exclude overflows by default (used in projections)
 159   if(fExclOffEntriesInProj){
 160     for(Int_t i=0;i<fNVar;i++){
 161       if(i!=ivar1 && i!=ivar2)
 162         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
 163     }
 164   }
 165   TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
 166   return hist;
 167
 168 }
 169 //___________________________________________________________________
 170 TH3D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
 171 {
 172   //
 173   // Make a 3D projection along variables ivar1 & ivar2 & ivar3
 174   //
 175   //exclude overflows by default (used in projections)
 176   if(fExclOffEntriesInProj){
 177     for(Int_t i=0;i<fNVar;i++){
 178       if(i!=ivar1 && i!=ivar2 && i!=ivar3)
 179         fData->GetAxis(i)->SetBit(TAxis::kAxisRange);
 180     }
 181   }
 182
 183   TH3D *hist=fData->Projection(ivar1,ivar2,ivar3);
 184   return hist;
 185
 186 }
 187
 188 //____________________________________________________________________
 189 Float_t AliCFGridSparse::GetOverFlows(Int_t ivar) const
 190 {
 191   //
 192   // Returns exclusive overflows in variable ivar
 193   //
 194   Int_t* bin = new Int_t[fNVar];
 195   memset(bin, 0, sizeof(Int_t) * fNVar);
 196   Float_t ovfl=0.;
 197   for (Long64_t i = 0; i < fData->GetNbins(); ++i) {
 198     Double_t v = fData->GetBinContent(i, bin);
 199     Bool_t add=kTRUE;
 200     for(Int_t j=0;j<fNVar;j++){
 201       if(ivar==j)continue;
 202       if((bin[j]==0) || (bin[j]==fNVarBins[j]+1))add=kFALSE;
 203     }
 204     if(bin[ivar]==fNVarBins[ivar]+1 && add) ovfl+=v;
 205   }
 206
 207   delete[] bin;
 208   return ovfl;
 209 }
 210
 211 //____________________________________________________________________
 212 Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar) const
 213 {
 214   //
 215   // Returns exclusive overflows in variable ivar
 216   //
 217   Int_t* bin = new Int_t[fNVar];
 218   memset(bin, 0, sizeof(Int_t) * fNVar);
 219   Float_t unfl=0.;
 220   for (Long64_t i = 0; i < fData->GetNbins(); ++i) {
 221     Double_t v = fData->GetBinContent(i, bin);
 222     Bool_t add=kTRUE;
 223     for(Int_t j=0;j<fNVar;j++){
 224       if(ivar==j)continue;
 225       if((bin[j]==0) || (bin[j]==fNVarBins[j]+1))add=kFALSE;
 226     }
 227     if(bin[ivar]==0 && add) unfl+=v;
 228   }
 229
 230   delete[] bin;
 231   return unfl;
 232 }
 233
 234
 235 //____________________________________________________________________
 236 Float_t AliCFGridSparse::GetEntries() const
 237 {
 238   //
 239   // total entries (including overflows and underflows)
 240   //
 241
 242   return fData->GetEntries();
 243 }
 244
 245 //____________________________________________________________________
 246 Float_t AliCFGridSparse::GetElement(Int_t index) const
 247 {
 248   //
 249   // Returns content of grid element index according to the
 250   // linear indexing in AliCFFrame
 251   //
 252   Int_t *bin = new Int_t[fNVar];
 253   GetBinIndex(index, bin);
 254   for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1; //consistency with AliCFGrid
 255   Float_t val=GetElement(fIndex);
 256   delete [] bin;
 257   return val;
 258
 259 }
 260 //____________________________________________________________________
 261 Float_t AliCFGridSparse::GetElement(Int_t *bin) const
 262 {
 263   //
 264   // Get the content in a bin corresponding to a set of bin indexes
 265   //
 266   return fData->GetBinContent(bin);
 267
 268 }
 269 //____________________________________________________________________
 270 Float_t AliCFGridSparse::GetElement(Double_t *var) const
 271 {
 272   //
 273   // Get the content in a bin corresponding to a set of input variables
 274   //
 275
 276   Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empty)
 277   if(index<0){
 278     return 0.;
 279   }else{
 280     return fData->GetBinContent(index);
 281   }
 282 }
 283
 284 //____________________________________________________________________
 285 Float_t AliCFGridSparse::GetElementError(Int_t index) const
 286 {
 287   //
 288   // Returns the error on the content of a bin according to a linear
 289   // indexing in AliCFFrame
 290   //
 291
 292   Int_t *bin = new Int_t[fNVar];
 293   GetBinIndex(index, bin);
 294   for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1; //consistency with AliCFGrid
 295   Float_t val=GetElementError(fIndex);
 296   delete [] bin;
 297   return val;
 298
 299 }
 300 //____________________________________________________________________
 301 Float_t AliCFGridSparse::GetElementError(Int_t *bin) const
 302 {
 303  //
 304   // Get the error in a bin corresponding to a set of bin indexes
 305   //
 306   return fData->GetBinError(bin);
 307
 308 }
 309 //____________________________________________________________________
 310 Float_t AliCFGridSparse::GetElementError(Double_t *var) const
 311 {
 312   //
 313   // Get the error in a bin corresponding to a set of input variables
 314   //
 315
 316   Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empy)
 317   if(index<0){
 318     return 0.;
 319   }else{
 320     return fData->GetBinError(index);
 321   }
 322 }
 323
 324
 325 //____________________________________________________________________
 326 void AliCFGridSparse::SetElement(Int_t index, Float_t val)
 327 {
 328   //
 329   // Sets grid element iel to val (linear indexing) in AliCFFrame
 330   //
 331   Int_t *bin = new Int_t[fNVar];
 332   GetBinIndex(index, bin);
 333   for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1;
 334   SetElement(fIndex,val);
 335   delete [] bin;
 336 }
 337 //____________________________________________________________________
 338 void AliCFGridSparse::SetElement(Int_t *bin, Float_t val)
 339 {
 340   //
 341   // Sets grid element of bin indeces bin to val
 342   //
 343   fData->SetBinContent(bin,val);
 344 }
 345 //____________________________________________________________________
 346 void AliCFGridSparse::SetElement(Double_t *var, Float_t val)
 347 {
 348   //
 349   // Set the content in a bin to value val corresponding to a set of input variables
 350   //
 351   Long_t index=fData->GetBin(var); //THnSparse index: allocate the cell
 352   Int_t *bin = new Int_t[fNVar];
 353   fData->GetBinContent(index,bin); //trick to access the array of bins
 354   fData->SetBinContent(bin,val);
 355   delete [] bin;
 356
 357 }
 358
 359 //____________________________________________________________________
 360 void AliCFGridSparse::SetElementError(Int_t index, Float_t val)
 361 {
 362   //
 363   // Sets grid element iel error to val (linear indexing) in AliCFFrame
 364   //
 365   Int_t *bin = new Int_t[fNVar];
 366   GetBinIndex(index, bin);
 367   for(Int_t i=0;i<fNVar;i++)fIndex[i]=bin[i]+1;
 368   SetElementError(fIndex,val);
 369   delete [] bin;
 370 }
 371 //____________________________________________________________________
 372 void AliCFGridSparse::SetElementError(Int_t *bin, Float_t val)
 373 {
 374   //
 375   // Sets grid element error of bin indeces bin to val
 376   //
 377   fData->SetBinError(bin,val);
 378 }
 379 //____________________________________________________________________
 380 void AliCFGridSparse::SetElementError(Double_t *var, Float_t val)
 381 {
 382   //
 383   // Set the error in a bin to value val corresponding to a set of input variables
 384   //
 385   Long_t index=fData->GetBin(var); //THnSparse index
 386   Int_t *bin = new Int_t[fNVar];
 387   fData->GetBinContent(index,bin); //trick to access the array of bins
 388   fData->SetBinError(bin,val);
 389   delete [] bin;
 390 }
 391
 392 //____________________________________________________________________
 393 void AliCFGridSparse::SumW2()
 394 {
 395   //
 396   //set calculation of the squared sum of the weighted entries
 397   //
 398   if(!fSumW2){
 399     fData->CalculateErrors(kTRUE);
 400   }
 401
 402   fSumW2=kTRUE;
 403 }
 404
 405 //____________________________________________________________________
 406 void AliCFGridSparse::Add(AliCFVGrid* aGrid, Double_t c)
 407 {
 408   //
 409   //add aGrid to the current one
 410   //
 411
 412
 413   if(aGrid->GetNVar()!=fNVar){
 414     AliInfo("Different number of variables, cannot add the grids");
 415     return;
 416   }
 417   if(aGrid->GetNDim()!=fNDim){
 418     AliInfo("Different number of dimensions, cannot add the grids!");
 419     return;
 420   }
 421
 422   if(!fSumW2  && aGrid->GetSumW2())SumW2();
 423
 424
 425   fData->Add(((AliCFGridSparse*)aGrid)->GetGrid(),c);
 426
 427 }
 428
 429 //____________________________________________________________________
 430 void AliCFGridSparse::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
 431 {
 432   //
 433   //Add aGrid1 and aGrid2 and deposit the result into the current one
 434   //
 435
 436   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
 437     AliInfo("Different number of variables, cannot add the grids");
 438     return;
 439   }
 440   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
 441     AliInfo("Different number of dimensions, cannot add the grids!");
 442     return;
 443   }
 444
 445   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
 446
 447
 448   fData->Reset();
 449   fData->Add(((AliCFGridSparse*)aGrid1)->GetGrid(),c1);
 450   fData->Add(((AliCFGridSparse*)aGrid2)->GetGrid(),c2);
 451
 452 }
 453
 454 //____________________________________________________________________
 455 void AliCFGridSparse::Multiply(AliCFVGrid* aGrid, Double_t c)
 456 {
 457   //
 458   // Multiply aGrid to the current one
 459   //
 460
 461
 462   if(aGrid->GetNVar()!=fNVar){
 463     AliInfo("Different number of variables, cannot multiply the grids");
 464     return;
 465   }
 466   if(aGrid->GetNDim()!=fNDim){
 467     AliInfo("Different number of dimensions, cannot multiply the grids!");
 468     return;
 469   }
 470
 471   if(!fSumW2  && aGrid->GetSumW2())SumW2();
 472
 473   THnSparse *h= ((AliCFGridSparse*)aGrid)->GetGrid();
 474
 475   fData->Multiply(h);
 476   fData->Scale(c);
 477
 478 }
 479
 480 //____________________________________________________________________
 481 void AliCFGridSparse::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
 482 {
 483   //
 484   //Multiply aGrid1 and aGrid2 and deposit the result into the current one
 485   //
 486
 487   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
 488     AliInfo("Different number of variables, cannot multiply the grids");
 489     return;
 490   }
 491   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
 492     AliInfo("Different number of dimensions, cannot multiply the grids!");
 493     return;
 494   }
 495
 496   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
 497
 498
 499   fData->Reset();
 500   THnSparse *h1= ((AliCFGridSparse*)aGrid1)->GetGrid();
 501   THnSparse *h2= ((AliCFGridSparse*)aGrid2)->GetGrid();
 502   h2->Multiply(h1);
 503   h2->Scale(c1*c2);
 504   fData->Add(h2);
 505 }
 506
 507
 508
 509 //____________________________________________________________________
 510 void AliCFGridSparse::Divide(AliCFVGrid* aGrid, Double_t c)
 511 {
 512   //
 513   // Divide aGrid to the current one
 514   //
 515
 516
 517   if(aGrid->GetNVar()!=fNVar){
 518     AliInfo("Different number of variables, cannot divide the grids");
 519     return;
 520   }
 521   if(aGrid->GetNDim()!=fNDim){
 522     AliInfo("Different number of dimensions, cannot divide the grids!");
 523     return;
 524   }
 525
 526   if(!fSumW2  && aGrid->GetSumW2())SumW2();
 527
 528   THnSparse *h= ((AliCFGridSparse*)aGrid)->GetGrid();
 529
 530   fData->Divide(h);
 531   fData->Scale(c);
 532
 533 }
 534
 535 //____________________________________________________________________
 536 void AliCFGridSparse::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
 537 {
 538   //
 539   //Divide aGrid1 and aGrid2 and deposit the result into the current one
 540   //bynomial errors are supported
 541   //
 542
 543   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
 544     AliInfo("Different number of variables, cannot divide the grids");
 545     return;
 546   }
 547   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
 548     AliInfo("Different number of dimensions, cannot divide the grids!");
 549     return;
 550   }
 551
 552   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
 553
 554
 555   THnSparse *h1= ((AliCFGridSparse*)aGrid1)->GetGrid();
 556   THnSparse *h2= ((AliCFGridSparse*)aGrid2)->GetGrid();
 557   fData->Divide(h1,h2,c1,c2,option);
 558 }
 559
 560
 561 //____________________________________________________________________
 562 void AliCFGridSparse::Rebin(const Int_t* group)
 563 {
 564   //
 565   // rebin the grid according to Rebin() as in THnSparse
 566   // Please notice that the original number of bins on
 567   // a given axis has to be divisible by the rebin group.
 568   //
 569
 570   for(Int_t i=0;i<fNVar;i++){
 571     if(group[i]!=1)AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
 572   }
 573
 574   THnSparse *rebinned =fData->Rebin(group);
 575   fData->Reset();
 576   fData = rebinned;
 577
 578   //redefine the needed stuff
 579
 580   Int_t ndimTot=1;
 581   Int_t nbinTot=0;
 582
 583   //number of bins in each dimension, auxiliary variables
 584
 585   for(Int_t ivar=0;ivar<fNVar;ivar++){
 586     Int_t nbins = fData->GetAxis(ivar)->GetNbins();
 587     fNVarBins[ivar]=nbins;
 588     ndimTot*=fNVarBins[ivar];
 589     nbinTot+=(fNVarBins[ivar]+1);
 590     Int_t offset=0;
 591     for(Int_t i =0;i<ivar;i++)offset+=(fNVarBins[i]+1);
 592     fOffset[ivar]=offset;
 593     Int_t prod=1;
 594     for(Int_t i=0;i<ivar;i++)prod*=fNVarBins[i];
 595     fProduct[ivar]=prod;
 596   }
 597
 598   fNDim=ndimTot;
 599
 600   //now the array of bin limits
 601
 602   delete fVarBinLimits;
 603   fNVarBinLimits=nbinTot;
 604   fVarBinLimits=new Double_t[fNVarBinLimits];
 605
 606   for(Int_t ivar=0;ivar<fNVar;ivar++){
 607     Double_t low = fData->GetAxis(ivar)->GetXmin();
 608     Double_t high = fData->GetAxis(ivar)->GetXmax();
 609     const TArrayD *xbins = fData->GetAxis(ivar)->GetXbins();
 610     if (xbins->fN == 0){
 611       for(Int_t ibin=0;ibin<=fNVarBins[ivar];ibin++){
 612         fVarBinLimits[ibin+fOffset[ivar]] = low + ibin*(high-low)/((Double_t) fNVarBins[ivar]);
 613       }
 614     }
 615     else{
 616
 617       for(Int_t ibin=0;ibin<=fNVarBins[ivar];ibin++) {
 618         fVarBinLimits[ibin+fOffset[ivar]] = xbins->At(ibin);
 619       }
 620     }
 621   }
 622
 623 }
 624 //____________________________________________________________________
 625 void AliCFGridSparse::Copy(TObject& c) const
 626 {
 627   //
 628   // copy function
 629   //
 630   AliCFGridSparse& target = (AliCFGridSparse &) c;
 631
 632   if(fData)target.fData = fData;
 633 }
 634