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