CORRFW/AliCFGrid.cxx

   1 /* $Id$ */
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16 //---------------------------------------------------------------------//
  17 //                                                                     //
  18 // AliCFGrid Class                                                     //
  19 // Class to accumulate data on an N-dimensional grid, to be used       //
  20 // as input to get corrections for Reconstruction & Trigger efficiency //
  21 // The class uses a one-dimensional array of floats to store the grid  //
  22 // --Author : S.Arcelli                                                //
  23 //   Still to be done:                                                 //
  24 // --Implement methods to merge cells                                  //
  25 // --Interpolate among bins in a range                                 //
  26 // This implementation will be aventually replaced  byAliCFGridSparse  //
  27 //---------------------------------------------------------------------//
  28 //
  29 //
  30 #include <AliLog.h>
  31 #include "AliCFGrid.h"
  32 #include "TMath.h"
  33 #include "TROOT.h"
  34 #include "TH1D.h"
  35 #include "TH2D.h"
  36 #include "TH3D.h"
  37
  38 //____________________________________________________________________
  39 ClassImp(AliCFGrid)
  40
  41 //____________________________________________________________________
  42 AliCFGrid::AliCFGrid() :
  43   AliCFVGrid(),
  44   fNunflTot(0),
  45   fNovflTot(0),
  46   fNentriesTot(0),
  47   fNunfl(0x0),
  48   fNovfl(0x0),
  49   fData(0x0),
  50   fErr2(0x0)
  51 {
  52   // default constructor
  53 }
  54 //____________________________________________________________________
  55 AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title) :
  56   AliCFVGrid(name,title),
  57   fNunflTot(0),
  58   fNovflTot(0),
  59   fNentriesTot(0),
  60   fNunfl(0x0),
  61   fNovfl(0x0),
  62   fData(0x0),
  63   fErr2(0x0)
  64 {
  65   // default constructor
  66 }
  67
  68 //____________________________________________________________________
  69 AliCFGrid::AliCFGrid(const Char_t* name, const Char_t* title, const Int_t nVarIn, const Int_t * nBinIn, const Double_t *binLimitsIn) :
  70   AliCFVGrid(name,title,nVarIn,nBinIn,binLimitsIn),
  71   fNunflTot(0),
  72   fNovflTot(0),
  73   fNentriesTot(0),
  74   fNunfl(0x0),
  75   fNovfl(0x0),
  76   fData(0x0),
  77   fErr2(0x0)
  78 {
  79   //
  80   // main constructor
  81   //
  82
  83   //The over/underflows
  84   fNunfl=new Float_t[fNVar];
  85   fNovfl= new Float_t[fNVar];
  86
  87
  88   //Initialization
  89   fNunflTot =0;
  90   fNovflTot =0;
  91   fNentriesTot =0;
  92   for(Int_t j=0;j<fNVar;j++){
  93     fNunfl[j] =0;
  94     fNovfl[j] =0;
  95   }
  96
  97
  98   // the grid
  99
 100   fData = new Float_t[fNDim]; //num
 101
 102   //Initialization
 103
 104 }
 105
 106 //____________________________________________________________________
 107 AliCFGrid::AliCFGrid(const AliCFGrid& c) :
 108   AliCFVGrid(c),
 109   fNunflTot(0),
 110   fNovflTot(0),
 111   fNentriesTot(0),
 112   fNunfl(0x0),
 113   fNovfl(0x0),
 114   fData(0x0),
 115   fErr2(0x0)
 116 {
 117   //
 118   // copy constructor
 119   //
 120   ((AliCFGrid &)c).Copy(*this);
 121 }
 122
 123 //____________________________________________________________________
 124 AliCFGrid::~AliCFGrid()
 125 {
 126   //
 127   // destructor
 128   //
 129   if(fNunfl)delete fNunfl;
 130   if(fNovfl)delete fNovfl;
 131   if(fData)delete fData;
 132   if(fSumW2)delete fErr2;
 133
 134 }
 135 //____________________________________________________________________
 136 AliCFGrid &AliCFGrid::operator=(const AliCFGrid &c)
 137 {
 138   //
 139   // assigment operator
 140   //
 141   if (this != &c)
 142     ((AliCFGrid &) c).Copy(*this);
 143   return *this;
 144 }
 145 //____________________________________________________________________
 146 Float_t AliCFGrid::GetElement(Int_t bin) const
 147 {
 148   //
 149   // Returns grid element bin
 150   //
 151   if(bin>=fNDim){
 152     AliInfo(Form(" element index outside the grid, return -1"));
 153     return -1.;
 154   }
 155   return fData[bin];
 156 }
 157 //____________________________________________________________________
 158 Float_t AliCFGrid::GetElement(Int_t *bin) const
 159 {
 160  //
 161   // Get the content in a bin corresponding to a set of bin indexes
 162   //
 163   //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
 164   for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
 165   Int_t ind =GetBinIndex(fIndex);
 166   return GetElement(ind);
 167 }
 168 //____________________________________________________________________
 169 Float_t AliCFGrid::GetElement(Double_t *var) const
 170 {
 171   //
 172   // Get the content in a bin corresponding to a set of input variables
 173   //
 174   Int_t unfl=0;
 175   Int_t ovfl=0;
 176   for(Int_t i=0;i<fNVar;i++){
 177     Int_t nbins=fNVarBins[i]+1;
 178     Float_t *bins=new Float_t[nbins];
 179     for(Int_t ibin =0;ibin<nbins;ibin++){
 180      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 181     }
 182
 183     fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
 184
 185     //underflows
 186
 187     if(var[i] < bins[0]){
 188       unfl=1;
 189     }
 190
 191     //overflows
 192
 193     if(var[i] > bins[nbins-1]){
 194       ovfl=1;
 195     }
 196     delete [] bins;
 197   }
 198
 199
 200   //move to the TH/THnSparse convention in N-dim bin numbering
 201   for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
 202
 203   if(!(ovfl==1 || unfl==1)){
 204     return GetElement(fIndex);
 205   }
 206   else{
 207     AliInfo(Form(" input variables outside the grid, return -1"));
 208     return -1.;
 209   }
 210 }
 211 //____________________________________________________________________
 212 Float_t AliCFGrid::GetElementError(Int_t iel) const
 213 {
 214   //
 215   // Return the error on grid element iel
 216   //
 217   if(iel>=fNDim){
 218     AliInfo(Form(" element index outside the grid, return -1"));
 219     return -1.;
 220   }
 221   if(fSumW2)return TMath::Sqrt(fErr2[iel]);
 222   return TMath::Sqrt(fData[iel]);
 223 }
 224 //____________________________________________________________________
 225 Float_t AliCFGrid::GetElementError(Int_t *bin) const
 226 {
 227   //
 228   // Get the error in a bin corresponding to a set of bin indeces
 229   //
 230   //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
 231   for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
 232   Int_t ind =GetBinIndex(fIndex);
 233   return GetElementError(ind);
 234
 235 }
 236 //____________________________________________________________________
 237 Float_t AliCFGrid::GetElementError(Double_t *var) const
 238 {
 239   //
 240   // Get the error in a bin corresponding to a set of input variables
 241   //
 242   Int_t unfl=0;
 243   Int_t ovfl=0;
 244   for(Int_t i=0;i<fNVar;i++){
 245     Int_t nbins=fNVarBins[i]+1;
 246     Float_t *bins=new Float_t[nbins];
 247     for(Int_t ibin =0;ibin<nbins;ibin++){
 248      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 249     }
 250
 251     fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
 252     //underflows
 253
 254     if(var[i] < bins[0]){
 255       unfl=1;
 256     }
 257
 258     //overflows
 259
 260     if(var[i] > bins[nbins-1]){
 261       ovfl=1;
 262     }
 263     delete [] bins;
 264   }
 265
 266   //move to the TH/THnSparse convention in N-dim bin numbering
 267   for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
 268   if(!(ovfl==1 || unfl==1)){
 269     return GetElementError(fIndex);
 270   }
 271   else{
 272     AliInfo(Form(" input variables outside the grid, return -1"));
 273     return -1.;
 274   }
 275 }
 276 //____________________________________________________________________
 277 void AliCFGrid::SetElement(Int_t iel, Float_t val)
 278 {
 279   //
 280   // Sets grid element iel to val
 281   //
 282   if(iel>=fNDim){
 283     AliInfo(Form(" element index outside the grid, no value set"));
 284   }else {
 285     fData[iel]=val;
 286   }
 287 }
 288 //____________________________________________________________________
 289 void AliCFGrid::SetElement(Int_t *bin, Float_t val)
 290 {
 291   //
 292   // Sets grid element of bin indeces bin to val
 293   //
 294   //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
 295   for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
 296   Int_t ind =GetBinIndex(fIndex);
 297   SetElement(ind,val);
 298 }
 299 //____________________________________________________________________
 300 void AliCFGrid::SetElement(Double_t *var, Float_t val)
 301 {
 302   //
 303   // Set the content in a bin to value val corresponding to a set of input variables
 304   //
 305   Int_t unfl=0;
 306   Int_t ovfl=0;
 307   for(Int_t i=0;i<fNVar;i++){
 308     Int_t nbins=fNVarBins[i]+1;
 309     Float_t *bins=new Float_t[nbins];
 310     for(Int_t ibin =0;ibin<nbins;ibin++){
 311      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 312     }
 313
 314     fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
 315     //underflows
 316
 317     if(var[i] < bins[0]){
 318       unfl=1;
 319     }
 320
 321     //overflows
 322
 323     if(var[i] > bins[nbins-1]){
 324       ovfl=1;
 325     }
 326     delete [] bins;
 327   }
 328
 329   //move to the TH/THnSparse convention in N-dim bin numbering
 330   for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
 331   if(!(ovfl==1 || unfl==1)){
 332     SetElement(fIndex,val);
 333   }
 334   else{
 335     AliInfo(Form(" input variables outside the grid, no value set"));
 336   }
 337 }
 338 //____________________________________________________________________
 339 void AliCFGrid::SetElementError(Int_t iel, Float_t val)
 340 {
 341   //
 342   // Set squared error on grid element iel to val*val
 343   //
 344   if(iel>=fNDim){
 345     AliInfo(Form(" element index outside the grid, no value set"));
 346     return;
 347   }
 348   if(!fErr2)SumW2();
 349   fErr2[iel]=val*val;
 350 }
 351 //____________________________________________________________________
 352 void AliCFGrid::SetElementError(Int_t *bin, Float_t val)
 353 {
 354   //
 355   // Set squared error to val on grid element of bin indeces bin
 356   //
 357   //-1 is to move from TH/ThnSparse N-dim bin convention to one in AliCFFrame
 358   for(Int_t i=0;i<fNVar; i++)fIndex[i]=bin[i]-1;
 359   Int_t ind =GetBinIndex(fIndex);
 360   SetElementError(ind,val);
 361 }
 362 //____________________________________________________________________
 363 void AliCFGrid::SetElementError(Double_t *var, Float_t val)
 364 {
 365   //
 366   // Set squared error to val in a bin corresponding to a set of input variables
 367   //
 368   Int_t unfl=0;
 369   Int_t ovfl=0;
 370   for(Int_t i=0;i<fNVar;i++){
 371     Int_t nbins=fNVarBins[i]+1;
 372     Float_t *bins=new Float_t[nbins];
 373     for(Int_t ibin =0;ibin<nbins;ibin++){
 374      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 375     }
 376
 377     fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
 378     //underflows
 379
 380     if(var[i] < bins[0]){
 381       unfl=1;
 382     }
 383
 384     //overflows
 385
 386     if(var[i] > bins[nbins-1]){
 387       ovfl=1;
 388     }
 389     delete [] bins;
 390   }
 391
 392   //move to the TH/THnSparse convention in N-dim bin numbering
 393   for(Int_t i=0;i<fNVar; i++)fIndex[i]+=1;
 394
 395   if(!(ovfl==1 || unfl==1)){
 396     SetElementError(fIndex,val);
 397   }
 398   else{
 399     AliInfo(Form(" input variables outside the grid, no value set"));
 400   }
 401 }
 402 //____________________________________________________________________
 403 void AliCFGrid::Fill(Double_t *var, Double_t weight)
 404 {
 405
 406   //
 407   // Fill the grid,
 408   // given a set of values of the input variable,
 409   // with weight (by default w=1)
 410   //
 411
 412
 413   Int_t isunfl=0;
 414   Int_t isovfl=0;
 415   Int_t *unfl=new Int_t[fNVar];
 416   Int_t *ovfl=new Int_t[fNVar];
 417
 418   for(Int_t i=0;i<fNVar;i++){
 419     unfl[i]=0;
 420     ovfl[i]=0;
 421   }
 422
 423   for(Int_t i=0;i<fNVar;i++){
 424     Int_t nbins=fNVarBins[i]+1;
 425     Float_t *bins=new Float_t[nbins];
 426     for(Int_t ibin =0;ibin<nbins;ibin++){
 427      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 428     }
 429
 430     fIndex[i] = TMath::BinarySearch(nbins,bins,var[i]);
 431     //underflows
 432
 433     if(var[i] < bins[0]){
 434       unfl[i]=1;
 435       isunfl=1;
 436     }
 437
 438     //overflows
 439
 440     if(var[i] > bins[nbins-1]){
 441       ovfl[i]=1;
 442       isovfl=1;
 443     }
 444     delete [] bins;
 445   }
 446
 447   //exclusive under/overflows
 448
 449   for(Int_t i=0;i<fNVar;i++){
 450     Bool_t add=kTRUE;
 451     for(Int_t j=0;j<fNVar;j++){
 452       if(i==j)continue;
 453       if(!(unfl[j]==0 && ovfl[j]==0))add=kFALSE;
 454     }
 455     if(add && unfl[i]==1)fNunfl[i]++;
 456     if(add && ovfl[i]==1)fNovfl[i]++;
 457   }
 458
 459   delete [] unfl;
 460   delete [] ovfl;
 461
 462   // Total number of entries, overflows and underflows
 463
 464   fNentriesTot++;
 465   if(isunfl)fNunflTot++;
 466   if(isovfl)fNovflTot++;
 467
 468   //if not ovfl/unfl, fill the element
 469
 470   if(!(isovfl==1 || isunfl==1)){
 471     Int_t ind =GetBinIndex(fIndex);
 472     fData[ind]+=weight;
 473     if(fSumW2)fErr2[ind]+=(weight*weight);
 474   }
 475 }
 476 //____________________________________________________________________
 477 Float_t AliCFGrid::GetOverFlows( Int_t ivar) const {
 478   //
 479   // Get overflows in variable var
 480   //
 481   return fNovfl[ivar];
 482 }
 483 //____________________________________________________________________
 484 Float_t AliCFGrid::GetOverFlows() const {
 485   //
 486   // Get overflows
 487   //
 488   return fNovflTot;
 489 }
 490 //____________________________________________________________________
 491 Float_t AliCFGrid::GetUnderFlows( Int_t ivar) const {
 492   //
 493   // Get overflows in variable var
 494   //
 495   return fNunfl[ivar];
 496 }
 497 //____________________________________________________________________
 498 Float_t AliCFGrid::GetUnderFlows() const {
 499   //
 500   // Get overflows
 501   //
 502   return fNunflTot;
 503 }
 504 //____________________________________________________________________
 505 Float_t AliCFGrid::GetEntries( ) const {
 506   //
 507   // Get total entries (in grid + over/underflows)
 508   //
 509   return fNentriesTot;
 510 }
 511 //___________________________________________________________________
 512 TH1D *AliCFGrid::Project(Int_t ivar) const
 513 {
 514   //
 515   // Make a 1D projection along variable ivar
 516
 517
 518   Int_t nbins =fNVarBins[ivar];
 519   Float_t *bins = new Float_t[nbins+1];
 520   for (Int_t i=0;i<=fNVar;i++){
 521   }
 522   for(Int_t ibin =0;ibin<nbins+1;ibin++){
 523     bins[ibin] = fVarBinLimits[ibin+fOffset[ivar]];
 524   }
 525
 526   char pname[40];
 527   sprintf(pname,"%s%s_%i",GetName(),"_proj1D_var", ivar);
 528   char htitle[40];
 529   sprintf(htitle,"%s%s_%i",GetName(),"_proj1D_var", ivar);
 530
 531   TH1D *proj1D=0;
 532
 533   //check if a projection with identical name exist
 534   TObject *obj = gROOT->FindObject(pname);
 535   if (obj && obj->InheritsFrom("TH1D")) {
 536     proj1D = (TH1D*)obj;
 537     proj1D->Reset();
 538   }
 539
 540   if(!proj1D){
 541     proj1D =new TH1D(pname,htitle, nbins, bins);
 542   }
 543
 544   delete [] bins;
 545   Float_t sum=0;
 546   Float_t *data= new Float_t[nbins];
 547   Float_t *err= new Float_t[nbins];
 548
 549   for(Int_t ibin=0;ibin<nbins;ibin++)data[ibin]=0;
 550   for(Int_t ibin=0;ibin<nbins;ibin++)err[ibin]=0;
 551   for(Int_t iel=0;iel<fNDim;iel++){
 552       data[GetBinIndex(ivar,iel)]+=fData[iel];
 553       if(fSumW2)err[GetBinIndex(ivar,iel)]+=fErr2[iel];
 554   }
 555
 556   for(Int_t ibin =0;ibin<nbins;ibin++){
 557     proj1D->SetBinContent(ibin+1,data[ibin]);
 558     proj1D->SetBinError(ibin+1,TMath::Sqrt(data[ibin]));
 559     if(fSumW2)proj1D->SetBinError(ibin+1,TMath::Sqrt(err[ibin]));
 560     sum+=data[ibin];
 561   }
 562
 563   delete [] data;
 564   delete [] err;
 565   proj1D->SetBinContent(nbins+1,GetOverFlows(ivar));
 566   proj1D->SetBinContent(0,GetUnderFlows(ivar));
 567   proj1D->SetEntries(sum+GetUnderFlows(ivar)+GetOverFlows(ivar));
 568   return proj1D;
 569 }
 570
 571 //___________________________________________________________________
 572 TH2D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2) const
 573 {
 574   //
 575   // Make a 2D projection along variable ivar
 576
 577   Int_t nbins1 =fNVarBins[ivar1];
 578   Int_t nbins2 =fNVarBins[ivar2];
 579
 580   Float_t *bins1 = new Float_t[nbins1+1];
 581   Float_t *bins2 = new Float_t[nbins2+1];
 582
 583   for(Int_t ibin =0;ibin<nbins1+1;ibin++){
 584     bins1[ibin] = fVarBinLimits[ibin+fOffset[ivar1]];
 585   }
 586   for(Int_t ibin =0;ibin<nbins2+1;ibin++){
 587     bins2[ibin] = fVarBinLimits[ibin+fOffset[ivar2]];
 588   }
 589
 590   char pname[40];
 591   sprintf(pname,"%s%s_%i_%i",GetName(),"_proj2D_var",ivar1,ivar2);
 592   char htitle[40];
 593   sprintf(htitle,"%s%s_%i_%i",GetName(),"_proj2D_var",ivar1,ivar2);
 594
 595   TH2D *proj2D=0;
 596
 597   //check if a projection with identical name exist
 598   TObject *obj = gROOT->FindObject(pname);
 599   if (obj && obj->InheritsFrom("TH2D")) {
 600     proj2D = (TH2D*)obj;
 601     proj2D->Reset();
 602   }
 603
 604   if(!proj2D){
 605     proj2D =new TH2D(pname,htitle, nbins1, bins1,nbins2,bins2);
 606   }
 607
 608   delete [] bins1;
 609   delete [] bins2;
 610
 611
 612   Float_t sum=0;
 613   Float_t **data=new Float_t*[nbins1];
 614   Float_t *data2=new Float_t[nbins1*nbins2];
 615   Float_t **err=new Float_t*[nbins1];
 616   Float_t *err2=new Float_t[nbins1*nbins2];
 617   for(Int_t i=0;i<nbins1;i++)data[i] = data2+i*nbins2;
 618   for(Int_t i=0;i<nbins1;i++)err[i] = err2+i*nbins2;
 619
 620   for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
 621     for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
 622       data[ibin1][ibin2]=0;
 623       err[ibin1][ibin2]=0;
 624     }
 625   }
 626
 627   for(Int_t iel=0;iel<fNDim;iel++){
 628     data[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)]+=fData[iel];
 629     if(fSumW2)err[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)]+=fErr2[iel];
 630   }
 631
 632   for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
 633     for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
 634       proj2D->SetBinContent(ibin1+1,ibin2+1,data[ibin1][ibin2]);
 635       proj2D->SetBinError(ibin1+1,ibin2+1,TMath::Sqrt(data[ibin1][ibin2]));
 636       if(fSumW2)proj2D->SetBinError(ibin1+1,ibin2+1,TMath::Sqrt(err[ibin1][ibin2]));
 637       sum+=data[ibin1][ibin2];
 638     }
 639
 640   }
 641   delete data;
 642   delete data2;
 643   delete err;
 644   delete err2;
 645   proj2D->SetBinContent(0,nbins2/2,GetUnderFlows(ivar1));
 646   proj2D->SetBinContent(nbins1+1,nbins2/2,GetOverFlows(ivar1));
 647   proj2D->SetBinContent(nbins1/2,0,GetUnderFlows(ivar2));
 648   proj2D->SetBinContent(nbins1/2,nbins2+1,GetOverFlows(ivar2));
 649   proj2D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2));
 650   return proj2D;
 651 }
 652 //___________________________________________________________________
 653 TH3D *AliCFGrid::Project(Int_t ivar1, Int_t ivar2, Int_t ivar3) const
 654 {
 655   //
 656   // Make a 3D projection along variable ivar
 657
 658   Int_t nbins1 =fNVarBins[ivar1];
 659   Int_t nbins2 =fNVarBins[ivar2];
 660   Int_t nbins3 =fNVarBins[ivar3];
 661
 662   Float_t *bins1 = new Float_t[nbins1+1];
 663   Float_t *bins2 = new Float_t[nbins2+1];
 664   Float_t *bins3 = new Float_t[nbins3+1];
 665
 666   for(Int_t ibin =0;ibin<nbins1+1;ibin++){
 667     bins1[ibin] = fVarBinLimits[ibin+fOffset[ivar1]];
 668   }
 669   for(Int_t ibin =0;ibin<nbins2+1;ibin++){
 670     bins2[ibin] = fVarBinLimits[ibin+fOffset[ivar2]];
 671   }
 672   for(Int_t ibin =0;ibin<nbins3+1;ibin++){
 673     bins3[ibin] = fVarBinLimits[ibin+fOffset[ivar3]];
 674   }
 675
 676   char pname[40];
 677   sprintf(pname,"%s%s_%i_%i_%i",GetName(),"_proj3D_var",ivar1,ivar2,ivar3);
 678   char htitle[40];
 679   sprintf(htitle,"%s%s_%i_%i_%i",GetName(),"_proj3D_var",ivar1,ivar2,ivar3);
 680
 681   TH3D *proj3D=0;
 682
 683   //check if a projection with identical name exist
 684   TObject *obj = gROOT->FindObject(pname);
 685   if (obj && obj->InheritsFrom("TH3D")) {
 686     proj3D = (TH3D*)obj;
 687     proj3D->Reset();
 688   }
 689
 690   if(!proj3D){
 691     proj3D =new TH3D(pname,htitle, nbins1,bins1,nbins2,bins2,nbins3,bins3);
 692   }
 693
 694   delete [] bins1;
 695   delete [] bins2;
 696   delete [] bins3;
 697
 698
 699   Float_t sum=0;
 700   Float_t ***data=new Float_t**[nbins1];
 701   Float_t **data2=new Float_t*[nbins1*nbins2];
 702   Float_t *data3=new Float_t[nbins1*nbins2*nbins3];
 703   Float_t ***err=new Float_t**[nbins1];
 704   Float_t **err2=new Float_t*[nbins1*nbins2];
 705   Float_t *err3=new Float_t[nbins1*nbins2*nbins3];
 706   for(Int_t i=0;i<nbins1;i++)data[i] = data2+i*nbins2;
 707   for(Int_t i=0;i<nbins1;i++)err[i] = err2+i*nbins2;
 708   for(Int_t i=0;i<nbins1;i++){
 709     for(Int_t j=0;j<nbins2;j++){
 710       data[i][j] = data3+i*nbins2*nbins3+j*nbins3;
 711       err[i][j] = err3+i*nbins2*nbins3+j*nbins3;
 712     }
 713   }
 714   for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
 715     for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
 716       for(Int_t ibin3 =0;ibin3<nbins3;ibin3++){
 717         data[ibin1][ibin2][ibin3]=0;
 718         err[ibin1][ibin2][ibin3]=0;
 719       }
 720     }
 721   }
 722
 723   for(Int_t iel=0;iel<fNDim;iel++){
 724     data[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)][GetBinIndex(ivar3,iel)]+=fData[iel];
 725     if(fSumW2)err[GetBinIndex(ivar1,iel)][GetBinIndex(ivar2,iel)][GetBinIndex(ivar3,iel)]+=fErr2[iel];
 726   }
 727
 728   for(Int_t ibin1 =0;ibin1<nbins1;ibin1++){
 729     for(Int_t ibin2 =0;ibin2<nbins2;ibin2++){
 730       for(Int_t ibin3 =0;ibin3<nbins3;ibin3++){
 731         proj3D->SetBinContent(ibin1+1,ibin2+1,ibin3+1,data[ibin1][ibin2][ibin3]);
 732         proj3D->SetBinError(ibin1+1,ibin2+1,ibin3+1,TMath::Sqrt(data[ibin1][ibin2][ibin3]));
 733         if(fSumW2)proj3D->SetBinError(ibin1+1,ibin2+1,ibin3+1,TMath::Sqrt(err[ibin1][ibin2][ibin3]));
 734         sum+=data[ibin1][ibin2][ibin3];
 735       }
 736     }
 737   }
 738
 739   delete data;
 740   delete data2;
 741   delete data3;
 742   delete err;
 743   delete err2;
 744   delete err3;
 745   proj3D->SetEntries(sum+GetUnderFlows(ivar1)+GetOverFlows(ivar1)+GetUnderFlows(ivar2)+GetOverFlows(ivar2)+GetUnderFlows(ivar3)+GetOverFlows(ivar3));
 746   return proj3D;
 747 }
 748
 749 //___________________________________________________________________
 750 TH1D *AliCFGrid::Slice(Int_t ivar, Double_t *varMin, Double_t* varMax) const
 751 {
 752   //
 753   // Make a slice along variable ivar in range [varMin,varMax]
 754
 755
 756   Int_t nbins =fNVarBins[ivar];
 757   Float_t *bins = new Float_t[nbins+1];
 758   for (Int_t i=0;i<=fNVar;i++){
 759   }
 760   for(Int_t ibin =0;ibin<nbins+1;ibin++){
 761     bins[ibin] = fVarBinLimits[ibin+fOffset[ivar]];
 762   }
 763
 764   char pname[40];
 765   sprintf(pname,"%s%s_%i",GetName(),"_proj1D_var", ivar);
 766   char htitle[40];
 767   sprintf(htitle,"%s%s_%i",GetName(),"_proj1D_var", ivar);
 768
 769   TH1D *proj1D=0;
 770
 771   //check if a projection with identical name exist
 772   TObject *obj = gROOT->FindObject(pname);
 773   if (obj && obj->InheritsFrom("TH1D")) {
 774     proj1D = (TH1D*)obj;
 775     proj1D->Reset();
 776   }
 777
 778   if(!proj1D){
 779     proj1D =new TH1D(pname,htitle, nbins, bins);
 780   }
 781
 782   delete [] bins;
 783
 784
 785   Int_t *indexMin=new Int_t[fNVar];
 786   Int_t *indexMax=new Int_t[fNVar];
 787
 788
 789   //Find out the min and max bins
 790
 791   for(Int_t i=0;i<fNVar;i++){
 792     Float_t xmin=varMin[i]; // the min values
 793     Float_t xmax=varMax[i]; // the min values
 794     Int_t nbins=fNVarBins[i]+1;
 795     Float_t *bins=new Float_t[nbins];
 796     for(Int_t ibin =0;ibin<nbins;ibin++){
 797      bins[ibin] = fVarBinLimits[ibin+fOffset[i]];
 798     }
 799     indexMin[i] = TMath::BinarySearch(nbins,bins,xmin);
 800     indexMax[i] = TMath::BinarySearch(nbins,bins,xmax);
 801     delete [] bins;
 802   }
 803
 804   Float_t sum=0;
 805   Float_t *data= new Float_t[nbins];
 806   for(Int_t ibin=0;ibin<nbins;ibin++)data[ibin]=0;
 807
 808   Int_t *index= new Int_t[fNVar];
 809   Int_t ielmin=GetBinIndex(indexMin);
 810   Int_t ielmax=GetBinIndex(indexMax);
 811   for(Int_t iel=ielmin;iel<=ielmax;iel++){
 812     GetBinIndex(iel,index);
 813     Bool_t isIn=kTRUE;
 814     for (Int_t j=0;j<fNVar;j++){
 815       if(!(index[j]>=indexMin[j] && index[j]<=indexMax[j]))isIn=kFALSE;
 816       break;
 817     }
 818     if(isIn)data[GetBinIndex(ivar,iel)]+=fData[iel];
 819   }
 820
 821   delete [] index;
 822
 823
 824   for(Int_t ibin =0;ibin<nbins;ibin++){
 825     proj1D->SetBinContent(ibin+1,data[ibin]);
 826     proj1D->SetBinError(ibin+1,TMath::Sqrt(data[ibin]));
 827     sum+=data[ibin];
 828   }
 829
 830   delete [] data;
 831
 832   proj1D->SetEntries(sum);
 833   return proj1D;
 834 }
 835
 836
 837 //____________________________________________________________________
 838 void AliCFGrid::Add(AliCFVGrid* aGrid, Double_t c)
 839 {
 840   //
 841   //add aGrid to the current one
 842   //
 843
 844   if(aGrid->GetNVar()!=fNVar){
 845     AliInfo("Different number of variables, cannot add the grids");
 846     return;
 847   }
 848   if(aGrid->GetNDim()!=fNDim){
 849     AliInfo("Different number of dimensions, cannot add the grids!");
 850     return;
 851   }
 852
 853   if(!fSumW2  && aGrid->GetSumW2())SumW2();
 854
 855   for(Int_t iel=0;iel<fNDim;iel++){
 856     fData[iel]+=(c*aGrid->GetElement(iel));
 857     if(fSumW2){
 858       Float_t err=aGrid->GetElementError(iel);
 859       fErr2[iel]+=c*c*err*err;
 860     }
 861   }
 862
 863   //Add entries, overflows and underflows
 864
 865   fNentriesTot+= c*aGrid->GetEntries();
 866   fNunflTot+= c*aGrid->GetUnderFlows();
 867   fNovflTot+= c*aGrid->GetOverFlows();
 868   for(Int_t j=0;j<fNVar;j++){
 869     fNunfl[j]+= c*aGrid->GetUnderFlows(j);
 870     fNovfl[j]+= c*aGrid->GetUnderFlows(j);
 871   }
 872 }
 873 //____________________________________________________________________
 874 void AliCFGrid::Add(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2)
 875 {
 876   //
 877   //add aGrid1 and aGrid2
 878   //
 879
 880   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
 881     AliInfo("Different number of variables, cannot add the grids");
 882     return;
 883   }
 884   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
 885     AliInfo("Different number of dimensions, cannot add the grids!");
 886     return;
 887   }
 888
 889   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
 890
 891   Float_t cont1,cont2,err1,err2;
 892
 893   for(Int_t iel=0;iel<fNDim;iel++){
 894     cont1=aGrid1->GetElement(iel);
 895     cont2=aGrid2->GetElement(iel);
 896     SetElement(iel,c1*cont1+c2*cont2);
 897     if(fSumW2){
 898       err1=aGrid1->GetElementError(iel);
 899       err2=aGrid2->GetElementError(iel);
 900       SetElementError(iel,TMath::Sqrt(c1*c1*err1*err1+c2*c2*err2*err2));
 901     }
 902   }
 903
 904   //Add entries, overflows and underflows
 905
 906   fNentriesTot= c1*aGrid1->GetEntries()+c2*aGrid2->GetEntries();
 907   fNunflTot= c1*aGrid1->GetUnderFlows()+c2*aGrid2->GetUnderFlows();
 908   fNovflTot= c1*aGrid1->GetOverFlows()+c2*aGrid2->GetOverFlows();
 909   for(Int_t j=0;j<fNVar;j++){
 910     fNunfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
 911     fNovfl[j]= c1*aGrid1->GetUnderFlows(j)+c2*aGrid2->GetUnderFlows(j);
 912   }
 913 }
 914 //____________________________________________________________________
 915 void AliCFGrid::Multiply(AliCFVGrid* aGrid, Double_t c)
 916 {
 917   //
 918   //multiply grid aGrid by the current one
 919   //
 920
 921   if(aGrid->GetNVar()!=fNVar){
 922     AliInfo("Different number of variables, cannot multiply the grids");
 923     return;
 924   }
 925   if(aGrid->GetNDim()!=fNDim){
 926     AliInfo("Different number of dimensions, cannot multiply the grids!");
 927     return;
 928   }
 929
 930   if(!fSumW2  && aGrid->GetSumW2())SumW2();
 931
 932   Float_t cont1,cont2,err1,err2;
 933
 934   for(Int_t iel=0;iel<fNDim;iel++){
 935     cont1=GetElement(iel);
 936     cont2=c*aGrid->GetElement(iel);
 937     SetElement(iel,cont1*cont2);
 938     if(fSumW2){
 939       err1=GetElementError(iel);
 940       err2=aGrid->GetElementError(iel);
 941       SetElementError(iel,TMath::Sqrt(c*c*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
 942     }
 943   }
 944
 945   //Set entries to the number of bins, preserve original overflows and underflows
 946
 947   fNentriesTot=fNDim;
 948   fNunflTot=GetUnderFlows();
 949   fNovflTot=GetOverFlows();
 950   for(Int_t j=0;j<fNVar;j++){
 951     fNunfl[j]= GetUnderFlows(j);
 952     fNovfl[j]= GetUnderFlows(j);
 953   }
 954 }
 955 //____________________________________________________________________
 956 void AliCFGrid::Multiply(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1, Double_t c2)
 957 {
 958   //
 959   //multiply grids aGrid1 and aGrid2
 960   //
 961
 962   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
 963     AliInfo("Different number of variables, cannot multiply the grids");
 964     return;
 965   }
 966   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
 967     AliInfo("Different number of dimensions, cannot multiply the grids!");
 968     return;
 969   }
 970
 971   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
 972
 973   Float_t cont1,cont2,err1,err2;
 974   for(Int_t iel=0;iel<fNDim;iel++){
 975     cont1=c1*aGrid1->GetElement(iel);
 976     cont2=c2*aGrid2->GetElement(iel);
 977     SetElement(iel,cont1*cont2);
 978     if(fSumW2){
 979       err1=aGrid1->GetElementError(iel);
 980       err2=aGrid2->GetElementError(iel);
 981       SetElementError(iel,TMath::Sqrt(c1*c1*c2*c2*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)));
 982     }
 983   }
 984
 985   //Set entries to the number of bins, preserve original overflows and underflows
 986
 987   fNentriesTot=fNDim;
 988   fNunflTot=GetUnderFlows();
 989   fNovflTot=GetOverFlows();
 990   for(Int_t j=0;j<fNVar;j++){
 991     fNunfl[j]= GetUnderFlows(j);
 992     fNovfl[j]= GetUnderFlows(j);
 993   }
 994 }
 995 //____________________________________________________________________
 996 void AliCFGrid::Divide(AliCFVGrid* aGrid, Double_t c)
 997 {
 998   //
 999   //divide current grid by grid aGrid
1000   //
1001
1002   if(aGrid->GetNVar()!=fNVar){
1003     AliInfo("Different number of variables, cannot divide the grids");
1004     return;
1005   }
1006   if(aGrid->GetNDim()!=fNDim){
1007     AliInfo("Different number of dimensions, cannot divide the grids!");
1008     return;
1009   }
1010  if(!c){AliInfo(Form("c is %f, cannot divide!",c)); return;}
1011
1012   if(!fSumW2  && aGrid->GetSumW2())SumW2();
1013
1014   Float_t cont1,cont2,err1,err2,den;
1015   for(Int_t iel=0;iel<fNDim;iel++){
1016     cont1=GetElement(iel);
1017     cont2=aGrid->GetElement(iel);
1018     if(cont2)SetElement(iel,cont1/(c*cont2));
1019     else SetElement(iel,0);
1020     if(fSumW2){
1021       err1=GetElementError(iel);
1022       err2=aGrid->GetElementError(iel);
1023       if(!cont2){SetElementError(iel,0.); continue;}
1024       den=cont2*cont2*cont2*c*c;
1025       SetElementError(iel,TMath::Sqrt((cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
1026     }
1027   }
1028
1029   //Set entries to the number of bins, preserve original overflows and underflows
1030
1031   fNentriesTot=fNDim;
1032   fNunflTot=GetUnderFlows();
1033   fNovflTot=GetOverFlows();
1034   for(Int_t j=0;j<fNVar;j++){
1035     fNunfl[j]= GetUnderFlows(j);
1036     fNovfl[j]= GetUnderFlows(j);
1037   }
1038 }
1039 //____________________________________________________________________
1040 void AliCFGrid::Divide(AliCFVGrid* aGrid1, AliCFVGrid* aGrid2, Double_t c1,Double_t c2, Option_t *option)
1041 {
1042   //
1043   //divide grids aGrid1,aGrid2
1044   //
1045
1046   TString opt = option;
1047   opt.ToUpper();
1048
1049   if(fNVar!=aGrid1->GetNVar()|| fNVar!=aGrid2->GetNVar()){
1050     AliInfo("Different number of variables, cannot divide the grids");
1051     return;
1052   }
1053   if(fNDim!=aGrid1->GetNDim()|| fNDim!=aGrid2->GetNDim()){
1054     AliInfo("Different number of dimensions, cannot divide the grids!");
1055     return;
1056   }
1057   if(!c2){AliInfo(Form("c2 is %f, cannot divide!",c2)); return;}
1058
1059   if(!fSumW2  && (aGrid1->GetSumW2() || aGrid2->GetSumW2()))SumW2();
1060
1061   Float_t cont1,cont2,err1,err2,r,den;
1062
1063   for(Int_t iel=0;iel<fNDim;iel++){
1064     cont1=aGrid1->GetElement(iel);
1065     cont2=aGrid2->GetElement(iel);
1066     if(cont2)SetElement(iel,c1*cont1/(c2*cont2));
1067     else SetElement(iel,0);
1068      if(fSumW2){
1069       err1=aGrid1->GetElementError(iel);
1070       err2=aGrid2->GetElementError(iel);
1071       if(!cont2){SetElementError(iel,0.); continue;}
1072       if (opt.Contains("B")){
1073         if(cont1!=cont2){
1074           r=cont1/cont2;
1075           SetElementError(iel,TMath::Sqrt(TMath::Abs(((1.-2.*r)*err1*err1+r*r*err2*err2)/(cont2*cont2))));
1076         }else{
1077           SetElementError(iel,0.);
1078         }
1079       }else{
1080         den=cont2*cont2*cont2*cont2*c2*c2;
1081         SetElementError(iel,TMath::Sqrt(c1*c1*(cont2*cont2*err1*err1+cont1*cont1*err2*err2)/den));
1082       }
1083     }
1084   }
1085
1086   //Set entries to the number of bins, preserve original overflows and underflows
1087
1088   fNentriesTot=fNDim;
1089   fNunflTot=GetUnderFlows();
1090   fNovflTot=GetOverFlows();
1091   for(Int_t j=0;j<fNVar;j++){
1092     fNunfl[j]= GetUnderFlows(j);
1093     fNovfl[j]= GetUnderFlows(j);
1094   }
1095 }
1096 //____________________________________________________________________
1097 void AliCFGrid::SumW2()
1098 {
1099   //
1100   //set calculation of the squared sum of the weighted entries
1101   //
1102   if(!fSumW2){
1103     fErr2=new Float_t [fNDim];
1104     //init....
1105     for(Int_t iel=0;iel<fNDim;iel++){
1106       fErr2[iel]=fData[iel];
1107     }
1108   }
1109
1110   fSumW2=kTRUE;
1111 }
1112 //____________________________________________________________________
1113 void AliCFGrid::Copy(TObject& c) const
1114 {
1115   //
1116   // copy function
1117   //
1118   AliCFGrid& target = (AliCFGrid &) c;
1119
1120   target.fNunflTot = fNunflTot;
1121   target.fNovflTot = fNovflTot;
1122   target.fNentriesTot = fNentriesTot;
1123   if (fNunfl)
1124     target.fNunfl = fNunfl;
1125   if (fNovfl)
1126     target.fNovfl = fNovfl;
1127   if (fData)
1128     target.fData = fData;
1129   if (fErr2)
1130     target.fErr2 = fErr2;
1131
1132 }