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