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