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   AliCFFrame(),
  43   fSumW2(kFALSE),
  44   fData(0x0)
  45 {
  46   // default constructor
  47 }
  48 //____________________________________________________________________
  49 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
  50   AliCFFrame(name,title),
  51   fSumW2(kFALSE),
  52   fData(0x0)
  53 {
  54   // default constructor
  55 }
  56 //____________________________________________________________________
  57 AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, Int_t nVarIn, const Int_t * nBinIn) :
  58   AliCFFrame(name,title),
  59   fSumW2(kFALSE),
  60   fData(0x0)
  61 {
  62   //
  63   // main constructor
  64   //
  65
  66   fData = new THnSparseF(name,title,nVarIn,nBinIn);
  67 }
  68
  69 //____________________________________________________________________
  70 AliCFGridSparse::~AliCFGridSparse()
  71 {
  72   //
  73   // destructor
  74   //
  75   if (fData) delete fData;
  76 }
  77
  78 //____________________________________________________________________
  79 AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
  80   AliCFFrame(c),
  81   fSumW2(kFALSE),
  82   fData(0x0)
  83 {
  84   //
  85   // copy constructor
  86   //
  87   ((AliCFGridSparse &)c).Copy(*this);
  88 }
  89
  90 //____________________________________________________________________
  91 AliCFGridSparse& AliCFGridSparse::operator=(const AliCFGridSparse &c)
  92 {
  93   //
  94   // assigment operator
  95   //
  96   if (this != &c) c.Copy(*this);
  97   return *this;
  98 }
  99
 100 //____________________________________________________________________
 101 void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t min, Double_t max)
 102 {
 103   //
 104   // set a uniform binning for variable ivar
 105   //
 106   Int_t nBins = GetNBins(ivar);
 107   Double_t * array = new Double_t[nBins+1];
 108   for (Int_t iEdge=0; iEdge<=nBins; iEdge++) array[iEdge] = min + iEdge * (max-min)/nBins ;
 109   fData->SetBinEdges(ivar, array);
 110   delete [] array ;
 111 }
 112
 113 //____________________________________________________________________
 114 void AliCFGridSparse::SetBinLimits(Int_t ivar, const Double_t *array)
 115 {
 116   //
 117   // setting the arrays containing the bin limits
 118   //
 119   fData->SetBinEdges(ivar, array);
 120 }
 121
 122 //____________________________________________________________________
 123 void AliCFGridSparse::Fill(const Double_t *var, Double_t weight)
 124 {
 125   //
 126   // Fill the grid,
 127   // given a set of values of the input variable,
 128   // with weight (by default w=1)
 129   //
 130   fData->Fill(var,weight);
 131 }
 132
 133 //___________________________________________________________________
 134 TH1D *AliCFGridSparse::Project(Int_t ivar) const
 135 {
 136   //
 137   // Make a 1D projection along variable ivar
 138   //
 139
 140   TH1D *hist=fData->Projection(ivar);
 141   hist->SetXTitle(GetVarTitle(ivar));
 142   hist->SetName(Form("%s_proj-%s",GetName(),GetVarTitle(ivar)));
 143   hist->SetTitle(Form("%s: projection on %s",GetTitle(),GetVarTitle(ivar)));
 144   for (Int_t iBin=1; iBin<=GetNBins(ivar); iBin++) hist->GetXaxis()->SetBinLabel(iBin,GetAxis(ivar)->GetBinLabel(iBin));
 145   return hist;
 146 }
 147 //___________________________________________________________________
 148 TH2D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2) const
 149 {
 150   //
 151   // Make a 2D projection along variables ivar1 & ivar2
 152   //
 153
 154   TH2D *hist=fData->Projection(ivar2,ivar1); //notice inverted axis (THnSparse uses TH3 2d-projection convention...)
 155   hist->SetXTitle(GetVarTitle(ivar1));
 156   hist->SetYTitle(GetVarTitle(ivar2));
 157   hist->SetName(Form("%s_proj-%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
 158   hist->SetTitle(Form("%s: projection on %s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2)));
 159   for (Int_t iBin=1; iBin<=GetNBins(ivar1); iBin++) hist->GetXaxis()->SetBinLabel(iBin,GetAxis(ivar1)->GetBinLabel(iBin));
 160   for (Int_t iBin=1; iBin<=GetNBins(ivar2); iBin++) hist->GetYaxis()->SetBinLabel(iBin,GetAxis(ivar2)->GetBinLabel(iBin));
 161   return hist;
 162 }
 163 //___________________________________________________________________
 164 TH3D *AliCFGridSparse::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
 165 {
 166   //
 167   // Make a 3D projection along variables ivar1 & ivar2 & ivar3
 168   //
 169
 170   TH3D *hist=fData->Projection(ivar1,ivar2,ivar3);
 171   hist->SetXTitle(GetVarTitle(ivar1));
 172   hist->SetYTitle(GetVarTitle(ivar2));
 173   hist->SetZTitle(GetVarTitle(ivar3));
 174   hist->SetName(Form("%s_proj-%s,%s,%s",GetName(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
 175   hist->SetTitle(Form("%s: projection on %s-%s-%s",GetTitle(),GetVarTitle(ivar1),GetVarTitle(ivar2),GetVarTitle(ivar3)));
 176   for (Int_t iBin=1; iBin<=GetNBins(ivar1); iBin++) hist->GetXaxis()->SetBinLabel(iBin,GetAxis(ivar1)->GetBinLabel(iBin));
 177   for (Int_t iBin=1; iBin<=GetNBins(ivar2); iBin++) hist->GetYaxis()->SetBinLabel(iBin,GetAxis(ivar2)->GetBinLabel(iBin));
 178   for (Int_t iBin=1; iBin<=GetNBins(ivar3); iBin++) hist->GetZaxis()->SetBinLabel(iBin,GetAxis(ivar3)->GetBinLabel(iBin));
 179   return hist;
 180 }
 181
 182 //___________________________________________________________________
 183 AliCFGridSparse* AliCFGridSparse::Project(Int_t nVars, const Int_t* vars, const Double_t* varMin, const Double_t* varMax, Bool_t useBins) const
 184 {
 185   //
 186   // projects the grid on the nVars dimensions defined in vars.
 187   // axis ranges can be defined in arrays varMin, varMax
 188   // If useBins=true, varMin and varMax are taken as bin numbers
 189   //
 190
 191   // binning for new grid
 192   Int_t* bins = new Int_t[nVars];
 193   for (Int_t iVar=0; iVar<nVars; iVar++) {
 194     bins[iVar] = GetNBins(vars[iVar]);
 195   }
 196
 197   // create new grid sparse
 198   AliCFGridSparse* out = new AliCFGridSparse(fName,fTitle,nVars,bins);
 199
 200   //set the range in the THnSparse to project
 201   THnSparse* clone = ((THnSparse*)fData->Clone());
 202   if (varMin && varMax) {
 203     for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
 204       if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
 205       else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
 206     }
 207   }
 208   else AliInfo("Keeping same axis ranges");
 209
 210   out->SetGrid(clone->Projection(nVars,vars));
 211   delete [] bins;
 212   return out;
 213 }
 214
 215
 216 //____________________________________________________________________
 217 Float_t AliCFGridSparse::GetBinCenter(Int_t ivar, Int_t ibin) const
 218 {
 219   //
 220   // Returns the center of specified bin for variable axis ivar
 221   //
 222
 223   return (Float_t) fData->GetAxis(ivar)->GetBinCenter(ibin);
 224 }
 225
 226 //____________________________________________________________________
 227 Float_t AliCFGridSparse::GetBinSize(Int_t ivar, Int_t ibin) const
 228 {
 229   //
 230   // Returns the size of specified bin for variable axis ivar
 231   //
 232
 233   return (Float_t) fData->GetAxis(ivar)->GetBinUpEdge(ibin) - fData->GetAxis(ivar)->GetBinLowEdge(ibin);
 234 }
 235
 236 //____________________________________________________________________
 237 Float_t AliCFGridSparse::GetEntries() const
 238 {
 239   //
 240   // total entries (including overflows and underflows)
 241   //
 242
 243   return fData->GetEntries();
 244 }
 245
 246 //____________________________________________________________________
 247 Float_t AliCFGridSparse::GetElement(Long_t index) const
 248 {
 249   //
 250   // Returns content of grid element index
 251   //
 252
 253   return fData->GetBinContent(index);
 254 }
 255 //____________________________________________________________________
 256 Float_t AliCFGridSparse::GetElement(const Int_t *bin) const
 257 {
 258   //
 259   // Get the content in a bin corresponding to a set of bin indexes
 260   //
 261   return fData->GetBinContent(bin);
 262
 263 }
 264 //____________________________________________________________________
 265 Float_t AliCFGridSparse::GetElement(const Double_t *var) const
 266 {
 267   //
 268   // Get the content in a bin corresponding to a set of input variables
 269   //
 270
 271   Long_t index = fData->GetBin(var,kFALSE);
 272   if (index<0) return 0.;
 273   return fData->GetBinContent(index);
 274 }
 275
 276 //____________________________________________________________________
 277 Float_t AliCFGridSparse::GetElementError(Long_t index) const
 278 {
 279   //
 280   // Returns the error on the content
 281   //
 282
 283   return fData->GetBinContent(index);
 284 }
 285 //____________________________________________________________________
 286 Float_t AliCFGridSparse::GetElementError(const Int_t *bin) const
 287 {
 288  //
 289   // Get the error in a bin corresponding to a set of bin indexes
 290   //
 291   return fData->GetBinError(bin);
 292
 293 }
 294 //____________________________________________________________________
 295 Float_t AliCFGridSparse::GetElementError(const Double_t *var) const
 296 {
 297   //
 298   // Get the error in a bin corresponding to a set of input variables
 299   //
 300
 301   Long_t index=fData->GetBin(var,kFALSE); //this is the THnSparse index (do not allocate new cells if content is empy)
 302   if (index<0) return 0.;
 303   return fData->GetBinError(index);
 304 }
 305
 306 //____________________________________________________________________
 307 void AliCFGridSparse::SetElement(Long_t index, Float_t val)
 308 {
 309   //
 310   // Sets grid element value
 311   //
 312   Int_t* bin = new Int_t[GetNVar()];
 313   fData->GetBinContent(index,bin); //affects the bin coordinates
 314   SetElement(bin,val);
 315   delete [] bin ;
 316 }
 317
 318 //____________________________________________________________________
 319 void AliCFGridSparse::SetElement(const Int_t *bin, Float_t val)
 320 {
 321   //
 322   // Sets grid element of bin indeces bin to val
 323   //
 324   fData->SetBinContent(bin,val);
 325 }
 326 //____________________________________________________________________
 327 void AliCFGridSparse::SetElement(const Double_t *var, Float_t val)
 328 {
 329   //
 330   // Set the content in a bin to value val corresponding to a set of input variables
 331   //
 332   Long_t index=fData->GetBin(var,kTRUE); //THnSparse index: allocate the cell
 333   Int_t *bin = new Int_t[GetNVar()];
 334   fData->GetBinContent(index,bin); //trick to access the array of bins
 335   SetElement(bin,val);
 336   delete [] bin;
 337 }
 338
 339 //____________________________________________________________________
 340 void AliCFGridSparse::SetElementError(Long_t index, Float_t val)
 341 {
 342   //
 343   // Sets grid element iel error to val (linear indexing) in AliCFFrame
 344   //
 345   Int_t *bin = new Int_t[GetNVar()];
 346   fData->GetBinContent(index,bin);
 347   SetElementError(bin,val);
 348   delete [] bin;
 349 }
 350
 351 //____________________________________________________________________
 352 void AliCFGridSparse::SetElementError(const Int_t *bin, Float_t val)
 353 {
 354   //
 355   // Sets grid element error of bin indeces bin to val
 356   //
 357   fData->SetBinError(bin,val);
 358 }
 359 //____________________________________________________________________
 360 void AliCFGridSparse::SetElementError(const Double_t *var, Float_t val)
 361 {
 362   //
 363   // Set the error in a bin to value val corresponding to a set of input variables
 364   //
 365   Long_t index=fData->GetBin(var); //THnSparse index
 366   Int_t *bin = new Int_t[GetNVar()];
 367   fData->GetBinContent(index,bin); //trick to access the array of bins
 368   SetElementError(bin,val);
 369   delete [] bin;
 370 }
 371
 372 //____________________________________________________________________
 373 void AliCFGridSparse::SumW2()
 374 {
 375   //
 376   //set calculation of the squared sum of the weighted entries
 377   //
 378   if(!fSumW2){
 379     fData->CalculateErrors(kTRUE);
 380   }
 381   fSumW2=kTRUE;
 382 }
 383
 384 //____________________________________________________________________
 385 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid, Double_t c)
 386 {
 387   //
 388   //add aGrid to the current one
 389   //
 390
 391   if (aGrid->GetNVar() != GetNVar()){
 392     AliError("Different number of variables, cannot add the grids");
 393     return;
 394   }
 395
 396   if (!fSumW2  && aGrid->GetSumW2()) SumW2();
 397   fData->Add(aGrid->GetGrid(),c);
 398 }
 399
 400 //____________________________________________________________________
 401 void AliCFGridSparse::Add(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
 402 {
 403   //
 404   //Add aGrid1 and aGrid2 and deposit the result into the current one
 405   //
 406
 407   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
 408     AliInfo("Different number of variables, cannot add the grids");
 409     return;
 410   }
 411
 412   if (!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
 413
 414   fData->Reset();
 415   fData->Add(aGrid1->GetGrid(),c1);
 416   fData->Add(aGrid2->GetGrid(),c2);
 417 }
 418
 419 //____________________________________________________________________
 420 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid, Double_t c)
 421 {
 422   //
 423   // Multiply aGrid to the current one
 424   //
 425
 426   if (aGrid->GetNVar() != GetNVar()) {
 427     AliError("Different number of variables, cannot multiply the grids");
 428     return;
 429   }
 430
 431   if(!fSumW2  && aGrid->GetSumW2()) SumW2();
 432   THnSparse *h = aGrid->GetGrid();
 433   fData->Multiply(h);
 434   fData->Scale(c);
 435 }
 436
 437 //____________________________________________________________________
 438 void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
 439 {
 440   //
 441   //Multiply aGrid1 and aGrid2 and deposit the result into the current one
 442   //
 443
 444   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
 445     AliError("Different number of variables, cannot multiply the grids");
 446     return;
 447   }
 448
 449   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
 450
 451   fData->Reset();
 452   THnSparse *h1 = aGrid1->GetGrid();
 453   THnSparse *h2 = aGrid2->GetGrid();
 454   h2->Multiply(h1);
 455   h2->Scale(c1*c2);
 456   fData->Add(h2);
 457 }
 458
 459 //____________________________________________________________________
 460 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid, Double_t c)
 461 {
 462   //
 463   // Divide aGrid to the current one
 464   //
 465
 466   if (aGrid->GetNVar() != GetNVar()) {
 467     AliError("Different number of variables, cannot divide the grids");
 468     return;
 469   }
 470
 471   if (!fSumW2  && aGrid->GetSumW2()) SumW2();
 472
 473   THnSparse *h1 = aGrid->GetGrid();
 474   THnSparse *h2 = (THnSparse*)fData->Clone();
 475   fData->Divide(h2,h1);
 476   fData->Scale(c);
 477 }
 478
 479 //____________________________________________________________________
 480 void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2, Option_t *option)
 481 {
 482   //
 483   //Divide aGrid1 and aGrid2 and deposit the result into the current one
 484   //binomial errors are supported
 485   //
 486
 487   if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
 488     AliError("Different number of variables, cannot divide the grids");
 489     return;
 490   }
 491
 492   if (!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
 493
 494   THnSparse *h1= aGrid1->GetGrid();
 495   THnSparse *h2= aGrid2->GetGrid();
 496   fData->Divide(h1,h2,c1,c2,option);
 497 }
 498
 499
 500 //____________________________________________________________________
 501 void AliCFGridSparse::Rebin(const Int_t* group)
 502 {
 503   //
 504   // rebin the grid according to Rebin() as in THnSparse
 505   // Please notice that the original number of bins on
 506   // a given axis has to be divisible by the rebin group.
 507   //
 508
 509   for(Int_t i=0;i<GetNVar();i++){
 510     if (group[i]!=1) AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
 511   }
 512
 513   THnSparse *rebinned =fData->Rebin(group);
 514   fData->Reset();
 515   fData = rebinned;
 516 }
 517 //____________________________________________________________________
 518 void AliCFGridSparse::Scale(Long_t index, const Double_t *fact)
 519 {
 520   //
 521   //scale content of a certain cell by (positive) fact (with error)
 522   //
 523
 524   if (GetElement(index)==0 || fact[0]==0) return;
 525
 526   Double_t in[2], out[2];
 527   in[0]=GetElement(index);
 528   in[1]=GetElementError(index);
 529   GetScaledValues(fact,in,out);
 530   SetElement(index,out[0]);
 531   if (fSumW2) SetElementError(index,out[1]);
 532 }
 533 //____________________________________________________________________
 534 void AliCFGridSparse::Scale(const Int_t *bin, const Double_t *fact)
 535 {
 536   //
 537   //scale content of a certain cell by (positive) fact (with error)
 538   //
 539   if(GetElement(bin)==0 || fact[0]==0)return;
 540
 541   Double_t in[2], out[2];
 542   in[0]=GetElement(bin);
 543   in[1]=GetElementError(bin);
 544   GetScaledValues(fact,in,out);
 545   SetElement(bin,out[0]);
 546   if(fSumW2)SetElementError(bin,out[1]);
 547
 548 }
 549 //____________________________________________________________________
 550 void AliCFGridSparse::Scale(const Double_t *var, const Double_t *fact)
 551 {
 552   //
 553   //scale content of a certain cell by (positive) fact (with error)
 554   //
 555   if(GetElement(var)==0 || fact[0]==0)return;
 556
 557   Double_t in[2], out[2];
 558   in[0]=GetElement(var);
 559   in[1]=GetElementError(var);
 560   GetScaledValues(fact,in,out);
 561   SetElement(var,out[0]);
 562   if(fSumW2)SetElementError(var,out[1]);
 563
 564 }
 565 //____________________________________________________________________
 566 void AliCFGridSparse::Scale(const Double_t* fact)
 567 {
 568   //
 569   //scale contents of the whole grid by fact
 570   //
 571
 572   for (Long_t iel=0; iel<GetNFilledBins(); iel++) {
 573     Scale(iel,fact);
 574   }
 575 }
 576 //____________________________________________________________________
 577 Long_t AliCFGridSparse::GetEmptyBins() const {
 578   //
 579   // Get empty bins
 580   //
 581
 582   return (GetNBinsTotal() - GetNFilledBins()) ;
 583 }
 584
 585 //_____________________________________________________________________
 586 // Int_t AliCFGridSparse::GetEmptyBins( Double_t *varMin, Double_t* varMax ) const
 587 // {
 588 //   //
 589 //   // Get empty bins in a range specified by varMin and varMax
 590 //   //
 591
 592 //   Int_t *indexMin = new Int_t[GetNVar()];
 593 //   Int_t *indexMax = new Int_t[GetNVar()];
 594
 595 //   //Find out the min and max bins
 596
 597 //   for (Int_t i=0;i<GetNVar();i++) {
 598 //     Double_t xmin=varMin[i]; // the min values
 599 //     Double_t xmax=varMax[i]; // the min values
 600 //     Int_t nbins = fData->GetAxis(i)->GetNbins()+1;
 601 //     Double_t *bins=new Double_t[nbins];
 602 //     for(Int_t ibin =0;ibin<nbins;ibin++){
 603 //      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 604 //     }
 605 //     indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
 606 //     indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
 607 //     delete [] bins;
 608 //   }
 609
 610 //   Int_t val=0;
 611 //   for(Int_t i=0;i<fNDim;i++){
 612 //     for (Int_t j=0;j<GetNVar();j++)fIndex[j]=GetBinIndex(j,i);
 613 //     Bool_t isIn=kTRUE;
 614 //     for (Int_t j=0;j<GetNVar();j++){
 615 //       if(!(fIndex[j]>=indexMin[j] && fIndex[j]<=indexMax[j]))isIn=kFALSE;
 616 //     }
 617 //     if(isIn && GetElement(i)<=0)val++;
 618 //   }
 619 //   AliInfo(Form(" the empty bins = %i ",val));
 620
 621 //   delete [] indexMin;
 622 //   delete [] indexMax;
 623 //   return val;
 624 // }
 625
 626 //____________________________________________________________________
 627 Int_t AliCFGridSparse::CheckStats(Double_t thr) const
 628 {
 629   //
 630   // Count the cells below a certain threshold
 631   //
 632   Int_t ncellsLow=0;
 633   for (Int_t i=0; i<GetNBinsTotal(); i++) {
 634     if (GetElement(i)<thr) ncellsLow++;
 635   }
 636   return ncellsLow;
 637 }
 638
 639 //_____________________________________________________________________
 640 Double_t AliCFGridSparse::GetIntegral() const
 641 {
 642   //
 643   // Get full Integral
 644   //
 645   return fData->ComputeIntegral();
 646 }
 647
 648 //_____________________________________________________________________
 649 // Double_t AliCFGridSparse::GetIntegral(Int_t *binMin, Int_t* binMax ) const
 650 // {
 651 //   //
 652 //   // Get Integral in a range of bin indeces (extremes included)
 653 //   //
 654
 655 //   Double_t val=0;
 656
 657 //   for(Int_t i=0;i<GetNVar();i++){
 658 //     if(binMin[i]<1)binMin[i]=1;
 659 //     if(binMax[i]>fNVarBins[i])binMax[i]=fNVarBins[i];
 660 //     if((binMin[i]>binMax[i])){
 661 //       AliInfo(Form(" Bin indices in variable %i in reverse order, please check!", i));
 662 //       return val;
 663 //     }
 664 //   }
 665 //   val=GetSum(0,binMin,binMax);
 666
 667 //   return val;
 668 // }
 669
 670 //_____________________________________________________________________
 671 // Double_t AliCFGridSparse::GetIntegral(Double_t *varMin, Double_t* varMax ) const
 672 // {
 673 //   //
 674 //   // Get Integral in a range (extremes included)
 675 //   //
 676
 677 //   Int_t *indexMin=new Int_t[GetNVar()];
 678 //   Int_t *indexMax=new Int_t[GetNVar()];
 679
 680 //   //Find out the min and max bins
 681
 682 //   for(Int_t i=0;i<GetNVar();i++){
 683 //     Double_t xmin=varMin[i]; // the min values
 684 //     Double_t xmax=varMax[i]; // the min values
 685 //     Int_t nbins=fNVarBins[i]+1;
 686 //     Double_t *bins=new Double_t[nbins];
 687 //     for(Int_t ibin =0;ibin<nbins;ibin++){
 688 //      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 689 //     }
 690 //     indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
 691 //     indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
 692 //     delete [] bins;
 693 //   }
 694
 695 //   //move to the TH/THnSparse convention in N-dim bin numbering
 696 //   for(Int_t i=0;i<GetNVar();i++){
 697 //     indexMin[i]+=1;
 698 //     indexMax[i]+=1;
 699 //   }
 700
 701 //   Double_t val=GetIntegral(indexMin,indexMax);
 702
 703 //   delete [] indexMin;
 704 //   delete [] indexMax;
 705
 706 //   return val;
 707 // }
 708
 709 // //_____________________________________________________________________
 710 // Double_t AliCFGridSparse::GetSum(Int_t ivar, Int_t *binMin, Int_t* binMax) const
 711 // {
 712 //   //
 713 //   // recursively add over nested loops....
 714 //   //
 715
 716 //   static Double_t val;
 717 //   if (ivar==0) val=0.;
 718 //   for(Int_t ibin=binMin[ivar]-1 ; ibin<=binMax[ivar]-1 ; ibin++){
 719 //     //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
 720 //     fIndex[ivar]=ibin;
 721 //     if(ivar<GetNVar()-1) {
 722 //       val=GetSum(ivar+1,binMin,binMax);
 723 //     }
 724 //     else {
 725 //       Int_t iel=GetBinIndex(fIndex);
 726 //       val+=GetElement(iel);
 727 //     }
 728 //   }
 729
 730 //   return val;
 731 // }
 732
 733 //____________________________________________________________________
 734 Long64_t AliCFGridSparse::Merge(TCollection* list)
 735 {
 736   //
 737   // Merge a list of AliCFGridSparse with this (needed for PROOF).
 738   // Returns the number of merged objects (including this).
 739   //
 740
 741   if (!list)
 742     return 0;
 743
 744   if (list->IsEmpty())
 745     return 1;
 746
 747   TIterator* iter = list->MakeIterator();
 748   TObject* obj;
 749
 750   Int_t count = 0;
 751   while ((obj = iter->Next())) {
 752     AliCFGridSparse* entry = dynamic_cast<AliCFGridSparse*> (obj);
 753     if (entry == 0)
 754       continue;
 755     this->Add(entry);
 756     count++;
 757   }
 758
 759   return count+1;
 760 }
 761
 762 //____________________________________________________________________
 763 void AliCFGridSparse::GetScaledValues(const Double_t *fact, const Double_t *in, Double_t *out) const{
 764   //
 765   // scale input *in and its error by (positive) fact (with error)
 766   // and erite it to *out
 767   //
 768   out[0]=in[0]*fact[0];
 769   out[1]=TMath::Sqrt(in[1]*in[1]/in[0]/in[0]
 770                      +fact[1]*fact[1]/fact[0]/fact[0])*out[0];
 771
 772 }
 773
 774 //____________________________________________________________________
 775 void AliCFGridSparse::Copy(TObject& c) const
 776 {
 777   //
 778   // copy function
 779   //
 780   AliCFFrame::Copy(c);
 781   AliCFGridSparse& target = (AliCFGridSparse &) c;
 782   target.fSumW2 = fSumW2 ;
 783   if (fData) {
 784     target.fData = (THnSparse*)fData->Clone();
 785   }
 786 }
 787
 788 //____________________________________________________________________
 789 TH1D* AliCFGridSparse::Slice(Int_t iVar, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
 790 {
 791   //
 792   // return a slice (1D-projection) on variable iVar while axis ranges are defined with varMin,varMax
 793   // arrays varMin and varMax contain the min and max values of each variable.
 794   // therefore varMin and varMax must have their dimensions equal to GetNVar()
 795   // If useBins=true, varMin and varMax are taken as bin numbers
 796
 797   THnSparse* clone = (THnSparse*)fData->Clone();
 798   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
 799     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
 800     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
 801   }
 802   TH1D* projection = clone->Projection(iVar);
 803   for (Int_t iBin=1; iBin<=GetNBins(iVar); iBin++) projection->GetXaxis()->SetBinLabel(iBin,GetAxis(iVar)->GetBinLabel(iBin));
 804   return projection ;
 805 }
 806
 807 //____________________________________________________________________
 808 TH2D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
 809 {
 810   //
 811   // return a slice (2D-projection) on variables iVar1 and iVar2 while axis ranges are defined with varMin,varMax
 812   // arrays varMin and varMax contain the min and max values of each variable.
 813   // therefore varMin and varMax must have their dimensions equal to GetNVar()
 814   // If useBins=true, varMin and varMax are taken as bin numbers
 815
 816   THnSparse* clone = (THnSparse*)fData->Clone();
 817   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
 818     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
 819     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
 820   }
 821   TH2D* projection = clone->Projection(iVar1,iVar2);
 822   for (Int_t iBin=1; iBin<=GetNBins(iVar1); iBin++) projection->GetXaxis()->SetBinLabel(iBin,GetAxis(iVar1)->GetBinLabel(iBin));
 823   for (Int_t iBin=1; iBin<=GetNBins(iVar2); iBin++) projection->GetYaxis()->SetBinLabel(iBin,GetAxis(iVar2)->GetBinLabel(iBin));
 824   return projection ;
 825 }
 826
 827 //____________________________________________________________________
 828 TH3D* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
 829 {
 830   //
 831   // return a slice (3D-projection) on variables iVar1, iVar2 and iVar3 while axis ranges are defined with varMin,varMax
 832   // arrays varMin and varMax contain the min and max values of each variable.
 833   // therefore varMin and varMax must have their dimensions equal to GetNVar()
 834   // If useBins=true, varMin and varMax are taken as bin numbers
 835
 836   THnSparse* clone = (THnSparse*)fData->Clone();
 837   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
 838     if (useBins)  clone->GetAxis(iAxis)->SetRange((Int_t)varMin[iAxis],(Int_t)varMax[iAxis]);
 839     else          clone->GetAxis(iAxis)->SetRangeUser(varMin[iAxis],varMax[iAxis]);
 840   }
 841   TH3D* projection = clone->Projection(iVar1,iVar2,iVar3);
 842   for (Int_t iBin=1; iBin<=GetNBins(iVar1); iBin++) projection->GetXaxis()->SetBinLabel(iBin,GetAxis(iVar1)->GetBinLabel(iBin));
 843   for (Int_t iBin=1; iBin<=GetNBins(iVar2); iBin++) projection->GetYaxis()->SetBinLabel(iBin,GetAxis(iVar2)->GetBinLabel(iBin));
 844   for (Int_t iBin=1; iBin<=GetNBins(iVar3); iBin++) projection->GetXaxis()->SetBinLabel(iBin,GetAxis(iVar3)->GetBinLabel(iBin));
 845   return projection ;
 846 }
 847
 848 //____________________________________________________________________
 849 void AliCFGridSparse::SetRangeUser(Int_t iVar, Double_t varMin, Double_t varMax) {
 850   //
 851   // set range of axis iVar.
 852   //
 853   fData->GetAxis(iVar)->SetRangeUser(varMin,varMax);
 854   AliWarning(Form("THnSparse axis %d range has been modified",iVar));
 855 }
 856
 857 //____________________________________________________________________
 858 void AliCFGridSparse::SetRangeUser(const Double_t *varMin, const Double_t *varMax) {
 859   //
 860   // set range of every axis. varMin and varMax must be of dimension GetNVar()
 861   //
 862   for (Int_t iAxis=0; iAxis<GetNVar() ; iAxis++) { // set new range for every axis
 863     SetRangeUser(iAxis,varMin[iAxis],varMax[iAxis]);
 864   }
 865   AliWarning("THnSparse axes ranges have been modified");
 866 }
 867
 868 //____________________________________________________________________
 869 void AliCFGridSparse::UseAxisRange(Bool_t b) const {
 870   for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) fData->GetAxis(iAxis)->SetBit(TAxis::kAxisRange,b);
 871 }
 872
 873 //____________________________________________________________________
 874 Float_t AliCFGridSparse::GetOverFlows(Int_t ivar, Bool_t exclusive) const
 875 {
 876   //
 877   // Returns overflows in variable ivar
 878   // Set 'exclusive' to true for an exclusive check on variable ivar
 879   //
 880   Int_t* bin = new Int_t[GetNVar()];
 881   memset(bin, 0, sizeof(Int_t) * GetNVar());
 882   Float_t ovfl=0.;
 883   for (Long64_t i = 0; i < fData->GetNbins(); i++) {
 884     Double_t v = fData->GetBinContent(i, bin);
 885     Bool_t add=kTRUE;
 886     if (exclusive) {
 887       for(Int_t j=0;j<GetNVar();j++){
 888         if(ivar==j)continue;
 889         if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
 890       }
 891     }
 892     if(bin[ivar]==GetNBins(ivar)+1 && add) ovfl+=v;
 893   }
 894
 895   delete[] bin;
 896   return ovfl;
 897 }
 898
 899 //____________________________________________________________________
 900 Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar, Bool_t exclusive) const
 901 {
 902   //
 903   // Returns exclusive overflows in variable ivar
 904   // Set 'exclusive' to true for an exclusive check on variable ivar
 905   //
 906   Int_t* bin = new Int_t[GetNVar()];
 907   memset(bin, 0, sizeof(Int_t) * GetNVar());
 908   Float_t unfl=0.;
 909   for (Long64_t i = 0; i < fData->GetNbins(); i++) {
 910     Double_t v = fData->GetBinContent(i, bin);
 911     Bool_t add=kTRUE;
 912     if (exclusive) {
 913       for(Int_t j=0;j<GetNVar();j++){
 914         if(ivar==j)continue;
 915         if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
 916       }
 917     }
 918     if(bin[ivar]==0 && add) unfl+=v;
 919   }
 920
 921   delete[] bin;
 922   return unfl;
 923 }
 924