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