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