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