PWG3/vertexingHF/AliMultiDimVector.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // Implementation of the class to store the number of signal     //
  21 // and background events in bins of the cut values               //
  22 // Origin:       Elena Bruna (bruna@to.infn.it)                  //
  23 // Updated:      Sergey Senyukov (senyukov@to.infn.it)           //
  24 //               Francesco Prino (prino@to.infn.it)              //
  25 // Last Updated: Giacomo Ortona (ortona@to.infn.it)              //
  26 //                                                               //
  27 ///////////////////////////////////////////////////////////////////
  28
  29 #include <fstream>
  30 #include <Riostream.h>
  31 #include "TH2.h"
  32 #include "AliMultiDimVector.h"
  33 #include "AliLog.h"
  34 #include "TString.h"
  35
  36 ClassImp(AliMultiDimVector)
  37 //___________________________________________________________________________
  38 AliMultiDimVector::AliMultiDimVector():TNamed("AliMultiDimVector","default"),
  39 fNVariables(0),
  40 fNPtBins(0),
  41 fVett(0),
  42 fNTotCells(0),
  43 fIsIntegrated(0)
  44 {
  45   // default constructor
  46 }
  47 //___________________________________________________________________________
  48 AliMultiDimVector::AliMultiDimVector(const char *name,const char *title, const Int_t nptbins, const Float_t* ptlimits, const Int_t npars,  const Int_t *nofcells,const Float_t *loosecuts, const Float_t *tightcuts, const TString *axisTitles):TNamed(name,title),
  49 fNVariables(npars),
  50 fNPtBins(nptbins),
  51 fVett(0),
  52 fNTotCells(0),
  53 fIsIntegrated(0){
  54 // standard constructor
  55   ULong64_t ntot=1;
  56   for(Int_t i=0;i<fNVariables;i++){
  57     ntot*=nofcells[i];
  58     fNCutSteps[i]=nofcells[i];
  59     if(loosecuts[i] == tightcuts[i]){
  60       if(loosecuts[i]!=0){
  61         printf("AliMultiDimVector::AliMultiDimVector: WARNING! same tight/loose variable for variable number %d. AliMultiDimVector with run with the following values: loose: %f; tight: %f\n",i,tightcuts[i]-0.1*tightcuts[i],tightcuts[i]);
  62         fMinLimits[i]=tightcuts[i]-0.1*tightcuts[i];
  63         fMaxLimits[i]=tightcuts[i];
  64       }else{
  65         fMinLimits[i]=0;
  66         fMaxLimits[i]=0.0001;
  67       }
  68         fGreaterThan[i]=kTRUE;
  69     }
  70     if(loosecuts[i] < tightcuts[i]){
  71       fMinLimits[i]=loosecuts[i];
  72       fMaxLimits[i]=tightcuts[i];
  73       fGreaterThan[i]=kTRUE;
  74     }else{
  75       fMinLimits[i]=tightcuts[i];
  76       fMaxLimits[i]=loosecuts[i];
  77       fGreaterThan[i]=kFALSE;
  78     }
  79     fAxisTitles[i]=axisTitles[i].Data();
  80   }
  81   fNTotCells=ntot*fNPtBins;
  82   fVett.Set(fNTotCells);
  83   for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=ptlimits[ipt];
  84   for(Int_t ipt=fNPtBins+1;ipt<fgkMaxNPtBins+1;ipt++) fPtLimits[ipt]=999.;
  85   for (ULong64_t j=0;j<fNTotCells;j++) fVett.AddAt(0,j);
  86 }
  87 //___________________________________________________________________________
  88 AliMultiDimVector::AliMultiDimVector(const AliMultiDimVector &mv):TNamed(mv.GetName(),mv.GetTitle()),
  89 fNVariables(mv.fNVariables),
  90 fNPtBins(mv.fNPtBins),
  91 fVett(0),
  92 fNTotCells(mv.fNTotCells),
  93 fIsIntegrated(mv.fIsIntegrated)
  94 {
  95   // copy constructor
  96   for(Int_t i=0;i<fNVariables;i++){
  97     fNCutSteps[i]=mv.GetNCutSteps(i);
  98     fMinLimits[i]=mv.GetMinLimit(i);
  99     fMaxLimits[i]=mv.GetMaxLimit(i);
 100     fGreaterThan[i]=mv.GetGreaterThan(i);
 101     fAxisTitles[i]=mv.GetAxisTitle(i);
 102   }
 103   fVett.Set(fNTotCells);
 104
 105   for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv.GetPtLimit(ipt);
 106   for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
 107 }
 108 //___________________________________________________________________________
 109 void AliMultiDimVector::CopyStructure(const AliMultiDimVector* mv){
 110 // Sets dimensions and limit from mv
 111   fNVariables=mv->GetNVariables();
 112   fNPtBins=mv->GetNPtBins();
 113   fNTotCells=mv->GetNTotCells();
 114   fIsIntegrated=mv->IsIntegrated();
 115   for(Int_t i=0;i<fNVariables;i++){
 116     fNCutSteps[i]=mv->GetNCutSteps(i);
 117     fMinLimits[i]=mv->GetMinLimit(i);
 118     fMaxLimits[i]=mv->GetMaxLimit(i);
 119     fGreaterThan[i]=mv->GetGreaterThan(i);
 120     fAxisTitles[i]=mv->GetAxisTitle(i);
 121   }
 122   for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv->GetPtLimit(ipt);
 123   fVett.Set(fNTotCells);
 124 }
 125 //______________________________________________________________________
 126 Bool_t AliMultiDimVector::GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const {
 127 // returns matrix element indices and Pt bin from global index
 128   if(globadd>=fNTotCells) return kFALSE;
 129   ULong64_t r=globadd;
 130   Int_t prod=1;
 131   Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
 132   for(Int_t k=0;k<fNVariables;k++) nOfCellsPlusLevel[k]=fNCutSteps[k];
 133   nOfCellsPlusLevel[fNVariables]=fNPtBins;
 134
 135   for(Int_t i=0;i<fNVariables+1;i++) prod*=nOfCellsPlusLevel[i];
 136   for(Int_t i=0;i<fNVariables+1;i++){
 137     prod/=nOfCellsPlusLevel[i];
 138     if(i<fNVariables) ind[i]=r/prod;
 139     else ptbin=r/prod;
 140     r=globadd%prod;
 141   }
 142   return kTRUE;
 143 }
 144 //______________________________________________________________________
 145 Bool_t AliMultiDimVector::GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const {
 146   Int_t ind[fgkMaxNVariables];
 147   Bool_t retcode=GetIndicesFromGlobalAddress(globadd,ind,ptbin);
 148   if(!retcode) return kFALSE;
 149   for(Int_t i=0;i<fNVariables;i++) cuts[i]=GetCutValue(i,ind[i]);
 150   return kTRUE;
 151 }
 152 //______________________________________________________________________
 153 ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(const Int_t *ind, Int_t ptbin) const {
 154   // Returns the global index of the cell in the matrix
 155   Int_t prod=1;
 156   ULong64_t elem=0;
 157   Int_t indexPlusLevel[fgkMaxNVariables+1];
 158   Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
 159   for(Int_t i=0;i<fNVariables;i++){
 160     indexPlusLevel[i]=ind[i];
 161     nOfCellsPlusLevel[i]=fNCutSteps[i];
 162   }
 163   indexPlusLevel[fNVariables]=ptbin;
 164   nOfCellsPlusLevel[fNVariables]=fNPtBins;
 165
 166   for(Int_t i=0;i<fNVariables+1;i++){
 167     prod=indexPlusLevel[i];
 168     if(i<fNVariables){
 169       for(Int_t j=i+1;j<fNVariables+1;j++){
 170         prod*=nOfCellsPlusLevel[j];
 171       }
 172     }
 173     elem+=prod;
 174   }
 175   return elem;
 176 }
 177 //______________________________________________________________________
 178 Bool_t AliMultiDimVector::GetIndicesFromValues(const Float_t *values, Int_t *ind) const {
 179   // Fills the array of matrix indices strating from variable values
 180   for(Int_t i=0;i<fNVariables;i++){
 181     if(fGreaterThan[i]){
 182       if(values[i]<GetMinLimit(i)) return kFALSE;
 183       ind[i]=(Int_t)((values[i]-fMinLimits[i])/GetCutStep(i));
 184       if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
 185     }else{
 186       if(values[i]>GetMaxLimit(i)) return kFALSE;
 187       ind[i]=(Int_t)((fMaxLimits[i]-values[i])/GetCutStep(i));
 188       if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
 189     }
 190   }
 191   return kTRUE;
 192 }
 193 //______________________________________________________________________
 194 ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(const Float_t *values, Int_t ptbin) const {
 195   // Returns the global index of the cell in the matrix
 196    Int_t ind[fgkMaxNVariables];
 197    Bool_t retcode=GetIndicesFromValues(values,ind);
 198    if(retcode) return GetGlobalAddressFromIndices(ind,ptbin);
 199    else{
 200      AliError("Values out of range");
 201      return fNTotCells+999;
 202    }
 203 }
 204 //_____________________________________________________________________________
 205 void AliMultiDimVector::MultiplyBy(Float_t factor){
 206   // multiply the AliMultiDimVector by a constant factor
 207   for(ULong64_t i=0;i<fNTotCells;i++){
 208     if(fVett.At(i)>0.)
 209       fVett.AddAt(fVett.At(i)*factor,i);
 210     else fVett.AddAt(-1,i);
 211   }
 212
 213 }
 214 //_____________________________________________________________________________
 215 void AliMultiDimVector::Multiply(const AliMultiDimVector* mv,Float_t factor){
 216   //  Sets AliMultiDimVector=mv*constant factor
 217   for(ULong64_t i=0;i<fNTotCells;i++){
 218     if(mv->GetElement(i)>0.)
 219       fVett.AddAt(mv->GetElement(i)*factor,i);
 220     else fVett.AddAt(-1,i);
 221   }
 222 }
 223 //_____________________________________________________________________________
 224 void AliMultiDimVector::Multiply(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
 225   //  Sets AliMultiDimVector=mv1*mv2
 226   for(ULong64_t i=0;i<fNTotCells;i++){
 227     if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
 228       fVett.AddAt(mv1->GetElement(i)*mv2->GetElement(i),i);
 229     else fVett.AddAt(-1,i);
 230   }
 231 }
 232 //_____________________________________________________________________________
 233 void AliMultiDimVector::Add(const AliMultiDimVector* mv){
 234   // Sums contents of mv to AliMultiDimVector
 235   if (mv->GetNTotCells()!=fNTotCells){
 236     AliError("Different dimension of the vectors!!");
 237   }else{
 238     for(ULong64_t i=0;i<fNTotCells;i++)
 239       if(mv->GetElement(i)>0. && fVett.At(i)>0.)
 240         fVett.AddAt(fVett.At(i)+mv->GetElement(i),i);
 241       else fVett.AddAt(-1,i);
 242   }
 243 }
 244 //_____________________________________________________________________________
 245 void AliMultiDimVector::Sum(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
 246   // Sets AliMultiDimVector=mv1+mv2
 247   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
 248     AliError("Different dimension of the vectors!!");
 249   }
 250   else{
 251     for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
 252       if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
 253         fVett.AddAt(mv1->GetElement(i)+mv2->GetElement(i),i);
 254       else fVett.AddAt(-1,i);
 255     }
 256   }
 257 }
 258 //_____________________________________________________________________________
 259 void AliMultiDimVector::LinearComb(const AliMultiDimVector* mv1, Float_t norm1, const AliMultiDimVector* mv2, Float_t norm2){
 260   // Sets AliMultiDimVector=n1*mv1+n2*mv2
 261   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
 262     AliError("Different dimension of the vectors!!");
 263   }
 264   else{
 265     for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
 266       if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
 267         fVett.AddAt(norm1*mv1->GetElement(i)+norm2*mv2->GetElement(i),i);
 268       else fVett.AddAt(-1,i);
 269     }
 270   }
 271 }
 272 //_____________________________________________________________________________
 273 void AliMultiDimVector::DivideBy(const AliMultiDimVector* mv){
 274   // Divide AliMulivector by mv
 275   if (mv->GetNTotCells()!=fNTotCells) {
 276     AliError("Different dimension of the vectors!!");
 277   }
 278   else{
 279     for(ULong64_t i=0;i<fNTotCells;i++)
 280       if(mv->GetElement(i)!=0 &&mv->GetElement(i)>0. && fVett.At(i)>0.)
 281         fVett.AddAt(fVett.At(i)/mv->GetElement(i),i);
 282       else fVett.AddAt(-1,i);
 283   }
 284
 285 }
 286 //_____________________________________________________________________________
 287 void AliMultiDimVector::Divide(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
 288   // Sets AliMultiDimVector=mv1/mv2
 289   if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
 290     AliError("Different dimension of the vectors!!");
 291   }
 292   else{
 293     for(ULong64_t i=0;i<mv1->GetNTotCells();i++)
 294       if(mv2->GetElement(i)!=0&& mv2->GetElement(i)>0.&& mv1->GetElement(i)>0.)
 295         {
 296           fVett.AddAt(mv1->GetElement(i)/mv2->GetElement(i),i);
 297         }
 298       else fVett.AddAt(-1,i);
 299   }
 300 }
 301 //_____________________________________________________________________________
 302 void AliMultiDimVector::Sqrt(){
 303   // Sqrt of elements of AliMultiDimVector
 304   for(ULong64_t i=0;i<fNTotCells;i++) {
 305     if(fVett.At(i)>=0) fVett.AddAt(TMath::Sqrt(fVett.At(i)),i);
 306     else {
 307       fVett.AddAt(-1,i);
 308     }
 309   }
 310 }
 311 //_____________________________________________________________________________
 312 void AliMultiDimVector::Sqrt(const AliMultiDimVector* mv){
 313   // Sets AliMultiDimVector=sqrt(mv)
 314   for(ULong64_t i=0;i<fNTotCells;i++)
 315     if(mv->GetElement(i)>=0) fVett.AddAt(TMath::Sqrt(mv->GetElement(i)),i);
 316     else fVett.AddAt(-1,i);
 317 }
 318 //_____________________________________________________________________________
 319 void AliMultiDimVector::FindMaximum(Float_t& maxValue, Int_t *ind , Int_t ptbin){
 320   // finds the element with maximum contents
 321   const ULong64_t nelem=fNTotCells/fNPtBins;
 322   TArrayF vett;
 323   vett.Set(nelem);
 324   ULong64_t runningAddress;
 325   for(ULong64_t i=0;i<nelem;i++){
 326     runningAddress=ptbin+i*fNPtBins;
 327     vett.AddAt(fVett[runningAddress],i);
 328   }
 329   maxValue=TMath::MaxElement(nelem,vett.GetArray());
 330   ULong64_t maxAddress=TMath::LocMax(nelem,vett.GetArray());
 331   ULong64_t maxGlobalAddress=ptbin+maxAddress*fNPtBins;
 332   Int_t checkedptbin;
 333   GetIndicesFromGlobalAddress(maxGlobalAddress,ind,checkedptbin);
 334 }
 335
 336 //_____________________________________________________________________________
 337 //Int_t* AliMultiDimVector::FindLocalMaximum(Float_t& maxValue, Bool_t *isFree,Int_t* indFixed, Int_t ptbin){
 338 Int_t* AliMultiDimVector::FindLocalMaximum(Float_t& maxValue, Int_t *numFixed,Int_t* indFixed, Int_t nfixed,Int_t ptbin){
 339   //return the elements with maximum content (maxValue) given fixed step for not free variables
 340   //numFixed[nfixed] is the indices of the fixed variables in the cuts array [fNVariables]={kTRUE,kTRUE,...,kFALSE,...,kFALSE,...,kTRUE}
 341   //indFixed[nfixed]={1,2} //nfixed is the number of false in isFree; indFixed contains the step for the i-th variable
 342   //!!take care of deleting the array of index returned!!
 343
 344   //  Int_t nfixed=0,nfree=0;
 345   //Int_t indtmp[fNVariables];
 346   if(nfixed>fNVariables)cout<<"AliMultiDimVector::FindLocalMaximum:ERROR! too many variables"<<endl;
 347   ULong64_t nelem=1;
 348   Int_t* indMax=new Int_t[fNVariables];
 349   //Get the number of fixed vars
 350   /*
 351   for (Int_t iv=0;iv<fNVariables;iv++){
 352     if(isFree[iv]){
 353       nfree++;
 354       nelem*=fNCutSteps[iv];
 355       indMax[iv]=0;
 356     }
 357     else {
 358       indMax[iv]=indFixed[nfixed];
 359       if(indFixed[nfixed]>=GetNCutSteps(iv)){
 360         indMax[iv]=0;
 361         cout<<"AliMultiDimVector::FindLocalMaximum:ERROR! called fixed ind "<< indFixed[nfixed]<<"  but "<<iv<<" var has only "<<GetNCutSteps(iv)<<" steps"<<endl;
 362       }
 363       nfixed++;
 364     }
 365   }
 366   */
 367   for (Int_t iv=0;iv<fNVariables;iv++)indMax[iv]=0;
 368   for(Int_t i=0;i<nfixed;i++){
 369     indMax[numFixed[i]]=indFixed[i];
 370   }
 371   //Get position of fixed vars
 372   /*
 373   Int_t fixedIndexes[nfixed];
 374   Int_t iforfixed=0;
 375   for (Int_t iv=0;iv<fNVariables;iv++){
 376     if(!isFree[iv]){
 377       fixedIndexes[iforfixed]=iv;
 378       iforfixed++;
 379     }
 380   }
 381   */
 382   TArrayF vett;
 383   vett.Set(nelem);
 384
 385   ULong64_t first=fNTotCells/fNPtBins*ptbin;
 386   ULong64_t last=first+fNTotCells/fNPtBins;
 387   Int_t dummyptbin;
 388
 389   maxValue=fVett[GetGlobalAddressFromIndices(indMax,ptbin)];
 390   Int_t tmpInd[fNVariables];
 391
 392   //loop on multidimvector global addresses
 393   for(ULong64_t iga=first;iga<last;iga++){
 394     GetIndicesFromGlobalAddress(iga,tmpInd,dummyptbin);
 395     Bool_t goodCell=kTRUE;
 396     for(Int_t ifix=0;ifix<nfixed&&goodCell;ifix++){
 397       //      if(indFixed[ifix]!=indMax[fixedIndexes[ifix]])goodCell=kFALSE;
 398       if(indFixed[ifix]!=tmpInd[numFixed[ifix]])goodCell=kFALSE;
 399     }
 400     if(goodCell){
 401       if(fVett[iga]>maxValue){
 402         maxValue=fVett[iga];
 403         //      GetIndicesFromGlobalAddress(iga,indMax,dummyptbin);
 404         for(Int_t inv=0;inv<fNVariables;inv++)indMax[inv]=tmpInd[inv];
 405       }
 406     }
 407   }
 408
 409   return indMax;
 410 }
 411
 412 //_____________________________________________________________________________
 413 TH2F*  AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, const Int_t* fixedVars, Int_t ptbin, Float_t norm){
 414   // Project the AliMultiDimVector on a 2D histogram
 415
 416   TString hisName=Form("hproj%s%dv%d",GetName(),secondVar,firstVar);
 417   TString hisTit=Form("%s vs. %s",fAxisTitles[secondVar].Data(),fAxisTitles[firstVar].Data());
 418   TH2F* h2=new TH2F(hisName.Data(),hisTit.Data(),fNCutSteps[firstVar],fMinLimits[firstVar],fMaxLimits[firstVar],fNCutSteps[secondVar],fMinLimits[secondVar],fMaxLimits[secondVar]);
 419
 420   Int_t index[fgkMaxNVariables];
 421   for(Int_t i=0;i<fNVariables;i++){
 422     index[i]=fixedVars[i];
 423   }
 424
 425   for(Int_t i=0;i<fNCutSteps[firstVar];i++){
 426     for(Int_t j=0;j<fNCutSteps[secondVar];j++){
 427       index[firstVar]=i;
 428       index[secondVar]=j;
 429       Float_t cont=GetElement(index,ptbin)/norm;
 430       Int_t bin1=i+1;
 431       if(!fGreaterThan[firstVar]) bin1=fNCutSteps[firstVar]-i;
 432       Int_t bin2=j+1;
 433       if(!fGreaterThan[secondVar]) bin2=fNCutSteps[secondVar]-j;
 434       h2->SetBinContent(bin1,bin2,cont);
 435     }
 436   }
 437   return h2;
 438 }
 439 //_____________________________________________________________________________
 440 void  AliMultiDimVector::GetIntegrationLimits(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
 441   // computes bin limits for integrating the AliMultiDimVector
 442   minbin=0;
 443   maxbin=0;
 444   if(iVar<fNVariables){
 445     minbin=iCell;
 446     maxbin=fNCutSteps[iVar]-1;
 447   }
 448 }
 449 //_____________________________________________________________________________
 450 void  AliMultiDimVector::GetFillRange(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
 451   // computes range of cells passing the cuts for FillAndIntegrate
 452   minbin=0;
 453   maxbin=0;
 454   if(iVar<fNVariables){
 455     minbin=0; // bin 0 corresponds to loose cuts
 456     maxbin=iCell;
 457   }
 458 }
 459 //_____________________________________________________________________________
 460 void AliMultiDimVector::Integrate(){
 461   // integrates the matrix
 462   if(fIsIntegrated){
 463     AliError("MultiDimVector already integrated");
 464     return;
 465   }
 466   TArrayF integral(fNTotCells);
 467   for(ULong64_t i=0;i<fNTotCells;i++) integral[i]=CountsAboveCell(i);
 468   for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]= integral[i];
 469   fIsIntegrated=kTRUE;
 470 }//_____________________________________________________________________________
 471 ULong64_t* AliMultiDimVector::GetGlobalAddressesAboveCuts(const Float_t *values, Int_t ptbin, Int_t& nVals) const{
 472   // fills an array with global addresses of cells passing the cuts
 473
 474   Int_t ind[fgkMaxNVariables];
 475   Bool_t retcode=GetIndicesFromValues(values,ind);
 476   if(!retcode){
 477     nVals=0;
 478     return 0x0;
 479   }
 480   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
 481   Int_t mink[fgkMaxNVariables];
 482   Int_t maxk[fgkMaxNVariables];
 483   Int_t size=1;
 484   for(Int_t i=0;i<fgkMaxNVariables;i++){
 485     GetFillRange(i,ind[i],mink[i],maxk[i]);
 486     size*=(maxk[i]-mink[i]+1);
 487   }
 488   ULong64_t* indexes=new ULong64_t[size];
 489   nVals=0;
 490   for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
 491     for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
 492       for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
 493         for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
 494           for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
 495             for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
 496               for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
 497                 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
 498                   for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
 499                     for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
 500                       Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
 501                       indexes[nVals++]=GetGlobalAddressFromIndices(currentBin,ptbin);
 502                     }
 503                   }
 504                 }
 505               }
 506             }
 507           }
 508         }
 509       }
 510     }
 511   }
 512   return indexes;
 513 }
 514 //_____________________________________________________________________________
 515 Float_t AliMultiDimVector::CountsAboveCell(ULong64_t globadd) const{
 516   // integrates the counts of cells above cell with address globadd
 517   Int_t ind[fgkMaxNVariables];
 518   Int_t ptbin;
 519   GetIndicesFromGlobalAddress(globadd,ind,ptbin);
 520   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
 521   Int_t mink[fgkMaxNVariables];
 522   Int_t maxk[fgkMaxNVariables];
 523   for(Int_t i=0;i<fgkMaxNVariables;i++){
 524     GetIntegrationLimits(i,ind[i],mink[i],maxk[i]);
 525   }
 526   Float_t sumcont=0.;
 527   for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
 528     for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
 529       for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
 530         for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
 531           for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
 532             for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
 533               for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
 534                 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
 535                   for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
 536                     for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
 537                       Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
 538                       sumcont+=GetElement(currentBin,ptbin);
 539                     }
 540                   }
 541                 }
 542               }
 543             }
 544           }
 545         }
 546       }
 547     }
 548   }
 549   return sumcont;
 550 }
 551 //_____________________________________________________________________________
 552 void AliMultiDimVector::Fill(Float_t* values, Int_t ptbin){
 553   // fills the cells of AliMultiDimVector corresponding to values
 554   if(fIsIntegrated){
 555     AliError("MultiDimVector already integrated -- Use FillAndIntegrate");
 556     return;
 557   }
 558   Int_t ind[fgkMaxNVariables];
 559   Bool_t retcode=GetIndicesFromValues(values,ind);
 560   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
 561   if(retcode) IncrementElement(ind,ptbin);
 562 }
 563 //_____________________________________________________________________________
 564 void AliMultiDimVector::FillAndIntegrate(Float_t* values, Int_t ptbin){
 565   // fills the cells of AliMultiDimVector passing the cuts
 566   // The number of nested loops must match fgkMaxNVariables!!!!
 567   fIsIntegrated=kTRUE;
 568   Int_t ind[fgkMaxNVariables];
 569   Bool_t retcode=GetIndicesFromValues(values,ind);
 570   if(!retcode) return;
 571   for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
 572   Int_t mink[fgkMaxNVariables];
 573   Int_t maxk[fgkMaxNVariables];
 574   for(Int_t i=0;i<fgkMaxNVariables;i++){
 575     GetFillRange(i,ind[i],mink[i],maxk[i]);
 576   }
 577   for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
 578     for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
 579       for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
 580         for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
 581           for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
 582             for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
 583               for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
 584                 for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
 585                   for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
 586                     for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
 587                       Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
 588                       IncrementElement(currentBin,ptbin);
 589                     }
 590                   }
 591                 }
 592               }
 593             }
 594           }
 595         }
 596       }
 597     }
 598   }
 599
 600 }
 601 //_____________________________________________________________________________
 602 void AliMultiDimVector::SuppressZeroBKGEffect(const AliMultiDimVector* mvBKG){
 603   // Sets to zero elements for which mvBKG=0
 604   for(ULong64_t i=0;i<fNTotCells;i++)
 605     if(mvBKG->GetElement(i)<0.00000001) fVett.AddAt(0,i);
 606 }
 607 //_____________________________________________________________________________
 608 AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBin){
 609   // sums the elements of pt bins between firstBin and lastBin
 610   if(firstBin<0 || lastBin>=fNPtBins || firstBin>=lastBin){
 611     AliError("Bad numbers of Pt bins to be shrinked");
 612     return 0;
 613   }
 614   Int_t nofcells[fgkMaxNVariables];
 615   Float_t loosecuts[fgkMaxNVariables];
 616   Float_t tightcuts[fgkMaxNVariables];
 617   TString axisTitles[fgkMaxNVariables];
 618   for(Int_t j=0;j<fgkMaxNVariables;j++) {
 619     nofcells[j]=0;
 620     loosecuts[j]=0.;
 621     tightcuts[j]=0.;
 622     axisTitles[j]="";
 623   }
 624   for(Int_t i=0;i<fNVariables;i++){
 625     nofcells[i]=fNCutSteps[i];
 626     if(fGreaterThan[i]){
 627       loosecuts[i]=fMinLimits[i];
 628       tightcuts[i]=fMaxLimits[i];
 629     }else{
 630       loosecuts[i]=fMaxLimits[i];
 631       tightcuts[i]=fMinLimits[i];
 632     }
 633     axisTitles[i]=fAxisTitles[i];
 634   }
 635   Int_t newNptbins=fNPtBins-(lastBin-firstBin);
 636   Float_t ptlimits[fgkMaxNPtBins+1];
 637   for(Int_t ipt=0; ipt<=firstBin;ipt++) ptlimits[ipt]=fPtLimits[ipt];
 638   for(Int_t ipt=firstBin+1; ipt<newNptbins+1;ipt++) ptlimits[ipt]=fPtLimits[ipt+(lastBin-firstBin)];
 639   AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),newNptbins,ptlimits,fNVariables,nofcells,loosecuts,tightcuts,axisTitles);
 640
 641   ULong64_t nOfPointsPerPtbin=fNTotCells/fNPtBins;
 642   ULong64_t addressOld,addressNew;
 643   Int_t npb,opb;
 644   for(npb=0;npb<firstBin;npb++){
 645     opb=npb;
 646     for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
 647       addressOld=opb+k*fNPtBins;
 648       addressNew=npb+k*newNptbins;
 649       shrinkedMV->SetElement(addressNew,fVett[addressOld]);
 650     }
 651   }
 652   npb=firstBin;
 653   for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
 654     Float_t summedValue=0.;
 655     for(opb=firstBin;opb<=lastBin;opb++){
 656       addressOld=opb+k*fNPtBins;
 657       summedValue+=fVett[addressOld];
 658     }
 659     addressNew=npb+k*newNptbins;
 660     shrinkedMV->SetElement(addressNew,summedValue);
 661   }
 662   for(npb=firstBin+1;npb<newNptbins;npb++){
 663     opb=npb+(lastBin-firstBin);
 664     for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
 665       addressOld=opb+k*fNPtBins;
 666       addressNew=npb+k*newNptbins;
 667       shrinkedMV->SetElement(addressNew,fVett[addressOld]);
 668     }
 669   }
 670   return shrinkedMV;
 671 }
 672 //_____________________________________________________________________________
 673 void AliMultiDimVector::SetNewLimits(Float_t* loose,Float_t* tight){
 674   for(Int_t i=0;i<fNVariables;i++){
 675     if(loose[i] < tight[i]){
 676       fMinLimits[i]=loose[i];
 677       fMaxLimits[i]=tight[i];
 678       fGreaterThan[i]=kTRUE;
 679     }else{
 680       fMinLimits[i]=tight[i];
 681       fMaxLimits[i]=loose[i];
 682       fGreaterThan[i]=kFALSE;
 683     }
 684   }
 685 }
 686 //_____________________________________________________________________________
 687 void AliMultiDimVector::SwapLimits(Int_t ivar){
 688   Float_t oldmin = fMinLimits[ivar];
 689   fMinLimits[ivar] = fMaxLimits[ivar];
 690   fMaxLimits[ivar] = oldmin;
 691   if(fGreaterThan[ivar])fGreaterThan[ivar]=kFALSE;
 692   else fGreaterThan[ivar]=kTRUE;
 693 }
 694 //_____________________________________________________________________________
 695 void AliMultiDimVector::PrintStatus(){
 696   //
 697   printf("Number of Pt bins       = %d\n",fNPtBins);
 698   printf("Limits of Pt bins       = ");
 699   for(Int_t ib=0;ib<fNPtBins+1;ib++) printf("%6.2f ",fPtLimits[ib]);
 700   printf("\n");
 701   printf("Number of cut variables = %d\n",fNVariables);
 702   for(Int_t iv=0;iv<fNVariables;iv++){
 703     printf("- Variable %d: %s\n",iv,fAxisTitles[iv].Data());
 704     printf("    Nsteps= %d Rage = %6.2f %6.2f\n",
 705            fNCutSteps[iv],fMinLimits[iv],fMaxLimits[iv]);
 706   }
 707 }