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