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