New classes and test macro for significance maximisation in multi-dimensional cuts...
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Feb 2009 17:13:40 +0000 (17:13 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Feb 2009 17:13:40 +0000 (17:13 +0000)
PWG3/vertexingHF/AliMultiDimVector.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliMultiDimVector.h [new file with mode: 0644]
PWG3/vertexingHF/AliSignificanceCalculator.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliSignificanceCalculator.h [new file with mode: 0644]
PWG3/vertexingHF/TestMultiVector.C [new file with mode: 0644]

diff --git a/PWG3/vertexingHF/AliMultiDimVector.cxx b/PWG3/vertexingHF/AliMultiDimVector.cxx
new file mode 100644 (file)
index 0000000..b6d61be
--- /dev/null
@@ -0,0 +1,527 @@
+/**************************************************************************
+ * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Implementation of the class to store the number of signal     //
+// and background events in bins of the cut values               //
+// Origin:       Elena Bruna (bruna@to.infn.it)                  //
+// Updated:      Sergey Senyukov (senyukov@to.infn.it)           //
+// Last updated: Francesco Prino (prino@to.infn.it)              //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+
+#include "TH2.h"
+#include "AliMultiDimVector.h"
+#include "AliLog.h"
+#include "TMath.h"
+#include "TString.h"
+
+ClassImp(AliMultiDimVector)
+//___________________________________________________________________________
+AliMultiDimVector::AliMultiDimVector():TNamed("AliMultiDimVector","default"),
+fNVariables(0),
+fNPtBins(0),
+fVett(0),
+fNTotCells(0),
+fIsIntegrated(0)
+{
+  // default constructor
+}
+//___________________________________________________________________________
+AliMultiDimVector::AliMultiDimVector(const char *name,const char *title, const Int_t npars, Int_t nptbins, Int_t *nofcells, Float_t *loosecuts, Float_t *tightcuts, TString *axisTitles):TNamed(name,title),
+fNVariables(npars),
+fNPtBins(nptbins),
+fVett(0),
+fNTotCells(0),
+fIsIntegrated(0){
+// standard constructor
+  ULong64_t ntot=1; 
+  for(Int_t i=0;i<fNVariables;i++){
+    ntot*=nofcells[i];
+    fNCutSteps[i]=nofcells[i];
+    if(loosecuts[i] <= tightcuts[i]){
+      fMinLimits[i]=loosecuts[i];
+      fMaxLimits[i]=tightcuts[i];
+      fGreaterThan[i]=kTRUE;
+    }else{
+      fMinLimits[i]=tightcuts[i];
+      fMaxLimits[i]=loosecuts[i];
+      fGreaterThan[i]=kFALSE;
+    }
+    fAxisTitles[i]=axisTitles[i].Data();
+  }
+  fNTotCells=ntot*fNPtBins;
+  fVett.Set(fNTotCells); 
+  for (ULong64_t j=0;j<fNTotCells;j++) fVett.AddAt(0,j);
+}
+//___________________________________________________________________________
+AliMultiDimVector::AliMultiDimVector(const AliMultiDimVector &mv):TNamed(mv.GetName(),mv.GetTitle()),
+fNVariables(mv.fNVariables),
+fNPtBins(mv.fNPtBins),
+fVett(0),
+fNTotCells(mv.fNTotCells),
+fIsIntegrated(mv.fIsIntegrated)
+{
+  // copy constructor
+  for(Int_t i=0;i<fNVariables;i++){
+    fNCutSteps[i]=mv.GetNCutSteps(i);
+    fMinLimits[i]=mv.GetMinLimit(i);
+    fMaxLimits[i]=mv.GetMaxLimit(i);
+    fGreaterThan[i]=mv.GetGreaterThan(i);
+    fAxisTitles[i]=mv.GetAxisTitle(i);
+  }
+  fVett.Set(fNTotCells); 
+  for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
+}
+//___________________________________________________________________________
+void AliMultiDimVector::CopyStructure(const AliMultiDimVector* mv){
+// Sets dimensions and limit from mv
+  fNVariables=mv->GetNVariables();
+  fNPtBins=mv->GetNPtBins();
+  fNTotCells=mv->GetNTotCells();
+  fIsIntegrated=mv->IsIntegrated();
+  for(Int_t i=0;i<fNVariables;i++){
+    fNCutSteps[i]=mv->GetNCutSteps(i);
+    fMinLimits[i]=mv->GetMinLimit(i);
+    fMaxLimits[i]=mv->GetMaxLimit(i);
+    fGreaterThan[i]=mv->GetGreaterThan(i);
+    fAxisTitles[i]=mv->GetAxisTitle(i);
+  }
+  fVett.Set(fNTotCells);
+  
+}
+//______________________________________________________________________
+Bool_t AliMultiDimVector::GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const {
+// returns matrix element indices and Pt bin from global index
+  if(globadd>=fNTotCells) return kFALSE;
+  ULong64_t r=globadd;
+  Int_t prod=1;
+  Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
+  for(Int_t k=0;k<fNVariables;k++) nOfCellsPlusLevel[k]=fNCutSteps[k];
+  nOfCellsPlusLevel[fNVariables]=fNPtBins;
+       
+  for(Int_t i=0;i<fNVariables+1;i++) prod*=nOfCellsPlusLevel[i];
+  for(Int_t i=0;i<fNVariables+1;i++){
+    prod/=nOfCellsPlusLevel[i];
+    if(i<fNVariables) ind[i]=r/prod;
+    else ptbin=r/prod;
+    r=globadd%prod;
+  }
+  return kTRUE;
+}
+//______________________________________________________________________
+Bool_t AliMultiDimVector::GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const {
+  Int_t ind[fgkMaxNVariables];
+  Bool_t retcode=GetIndicesFromGlobalAddress(globadd,ind,ptbin);
+  if(!retcode) return kFALSE;
+  for(Int_t i=0;i<fNVariables;i++) cuts[i]=GetCutValue(i,ind[i]);
+  return kTRUE;
+}
+//______________________________________________________________________
+ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(Int_t *ind, Int_t ptbin) const {
+  // Returns the global index of the cell in the matrix
+  Int_t prod=1;
+  ULong64_t elem=0;
+  Int_t indexPlusLevel[fgkMaxNVariables+1];
+  Int_t nOfCellsPlusLevel[fgkMaxNVariables+1];
+  for(Int_t i=0;i<fNVariables;i++){
+    indexPlusLevel[i]=ind[i];
+    nOfCellsPlusLevel[i]=fNCutSteps[i];
+  }
+  indexPlusLevel[fNVariables]=ptbin;
+  nOfCellsPlusLevel[fNVariables]=fNPtBins;
+       
+  for(Int_t i=0;i<fNVariables+1;i++){
+    prod=indexPlusLevel[i];
+    if(i<fNVariables){
+      for(Int_t j=i+1;j<fNVariables+1;j++){
+       prod*=nOfCellsPlusLevel[j];
+      }
+    }
+    elem+=prod;
+  }
+  return elem;
+}
+//______________________________________________________________________
+Bool_t AliMultiDimVector::GetIndicesFromValues(Float_t *values, Int_t *ind) const {
+  // Fills the array of matrix indices strating from variable values
+  for(Int_t i=0;i<fNVariables;i++){
+    if(fGreaterThan[i]){ 
+      if(values[i]<GetMinLimit(i)) return kFALSE;
+      ind[i]=(Int_t)((values[i]-fMinLimits[i])/GetCutStep(i));
+      if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
+    }else{
+      if(values[i]>GetMaxLimit(i)) return kFALSE;
+      ind[i]=(Int_t)((fMaxLimits[i]-values[i])/GetCutStep(i));
+      if(ind[i]>=GetNCutSteps(i)) ind[i]=GetNCutSteps(i)-1;
+    }
+  }
+  return kTRUE;
+}
+//______________________________________________________________________
+ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(Float_t *values, Int_t ptbin) const {
+  // Returns the global index of the cell in the matrix
+   Int_t ind[fgkMaxNVariables];
+   Bool_t retcode=GetIndicesFromValues(values,ind);
+   if(retcode) return GetGlobalAddressFromIndices(ind,ptbin);
+   else{
+     AliError("Values out of range");
+     return fNTotCells+999;
+   }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::MultiplyBy(Float_t factor){
+  // multiply the AliMultiDimVector by a constant factor
+  for(ULong64_t i=0;i<fNTotCells;i++){
+    if(fVett.At(i)!=-1)
+      fVett.AddAt(fVett.At(i)*factor,i);
+    else fVett.AddAt(-1,i);
+  }
+  
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Multiply(AliMultiDimVector* mv,Float_t factor){
+  //  Sets AliMultiDimVector=mv*constant factor
+  for(ULong64_t i=0;i<fNTotCells;i++){
+    if(mv->GetElement(i)!=-1)
+      fVett.AddAt(mv->GetElement(i)*factor,i);
+    else fVett.AddAt(-1,i); 
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Multiply(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+  //  Sets AliMultiDimVector=mv1*mv2
+  for(ULong64_t i=0;i<fNTotCells;i++){
+    if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
+      fVett.AddAt(mv1->GetElement(i)*mv2->GetElement(i),i);
+    else fVett.AddAt(-1,i); 
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Add(AliMultiDimVector* mv){
+  // Sums contents of mv to AliMultiDimVector
+  if (mv->GetNTotCells()!=fNTotCells){ 
+    AliError("Different dimension of the vectors!!");
+  }else{
+    for(ULong64_t i=0;i<fNTotCells;i++) 
+      if(mv->GetElement(i)!=-1 && fVett.At(i)!=-1)
+       fVett.AddAt(fVett.At(i)+mv->GetElement(i),i);
+      else fVett.AddAt(-1,i); 
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Sum(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+  // Sets AliMultiDimVector=mv1+mv2
+  if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
+    AliError("Different dimension of the vectors!!");
+  }
+  else{
+    for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
+      if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
+       fVett.AddAt(mv1->GetElement(i)+mv2->GetElement(i),i); 
+      else fVett.AddAt(-1,i); 
+    }
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::LinearComb(AliMultiDimVector* mv1, Float_t norm1, AliMultiDimVector* mv2, Float_t norm2){
+  // Sets AliMultiDimVector=n1*mv1+n2*mv2
+  if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
+    AliError("Different dimension of the vectors!!");
+  }
+  else{
+    for(ULong64_t i=0;i<mv1->GetNTotCells();i++) {
+      if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
+       fVett.AddAt(norm1*mv1->GetElement(i)+norm2*mv2->GetElement(i),i); 
+      else fVett.AddAt(-1,i); 
+    }
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::DivideBy(AliMultiDimVector* mv){
+  // Divide AliMulivector by mv
+  if (mv->GetNTotCells()!=fNTotCells) {
+    AliError("Different dimension of the vectors!!");
+  }
+  else{
+    for(ULong64_t i=0;i<fNTotCells;i++) 
+      if(mv->GetElement(i)!=0 &&mv->GetElement(i)!=-1 && fVett.At(i)!=-1)
+       fVett.AddAt(fVett.At(i)/mv->GetElement(i),i);
+      else fVett.AddAt(-1,i);
+  }
+
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Divide(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+  // Sets AliMultiDimVector=mv1/mv2
+  if (fNTotCells!=mv1->GetNTotCells()&&mv1->GetNTotCells()!=mv2->GetNTotCells()) {
+    AliError("Different dimension of the vectors!!");
+  }
+  else{
+    for(ULong64_t i=0;i<mv1->GetNTotCells();i++) 
+      if(mv2->GetElement(i)!=0&& mv2->GetElement(i)!=-1&& mv1->GetElement(i)!=-1)
+       {
+         fVett.AddAt(mv1->GetElement(i)/mv2->GetElement(i),i);
+       }
+      else fVett.AddAt(-1,i);
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Sqrt(){
+  // Sqrt of elements of AliMultiDimVector
+  for(ULong64_t i=0;i<fNTotCells;i++) {
+    if(fVett.At(i)>=0) fVett.AddAt(TMath::Sqrt(fVett.At(i)),i);
+    else {
+      fVett.AddAt(-1,i);
+    }
+  }
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::Sqrt(AliMultiDimVector* mv){
+  // Sets AliMultiDimVector=sqrt(mv)
+  for(ULong64_t i=0;i<fNTotCells;i++) 
+    if(mv->GetElement(i)>=0) fVett.AddAt(TMath::Sqrt(mv->GetElement(i)),i);
+    else fVett.AddAt(-1,i);
+}
+//_____________________________________________________________________________
+void AliMultiDimVector::FindMaximum(Float_t& maxValue, Int_t *ind , Int_t ptbin){
+  // finds the element with maximum contents
+  const ULong64_t nelem=fNTotCells/fNPtBins;
+  TArrayF vett;
+  vett.Set(nelem);
+  ULong64_t runningAddress;
+  for(ULong64_t i=0;i<nelem;i++){
+    runningAddress=ptbin+i*fNPtBins;
+    vett.AddAt(fVett[runningAddress],i);
+  }
+  maxValue=TMath::MaxElement(nelem,vett.GetArray());
+  ULong64_t maxAddress=TMath::LocMax(nelem,vett.GetArray());
+  ULong64_t maxGlobalAddress=ptbin+maxAddress*fNPtBins;
+  Int_t checkedptbin;
+  GetIndicesFromGlobalAddress(maxGlobalAddress,ind,checkedptbin);
+}
+
+//_____________________________________________________________________________
+TH2F*  AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, Int_t* fixedVars, Int_t ptbin, Float_t norm){
+  // Project the AliMultiDimVector on a 2D histogram
+
+  TString hisName=Form("hproj%s%dv%d",GetName(),secondVar,firstVar);
+  TString hisTit=Form("%s vs. %s",fAxisTitles[secondVar].Data(),fAxisTitles[firstVar].Data());
+  TH2F* h2=new TH2F(hisName.Data(),hisTit.Data(),fNCutSteps[firstVar],fMinLimits[firstVar],fMaxLimits[firstVar],fNCutSteps[secondVar],fMinLimits[secondVar],fMaxLimits[secondVar]);
+       
+  Int_t index[fgkMaxNVariables];
+  for(Int_t i=0;i<fNVariables;i++){
+    index[i]=fixedVars[i];
+  }
+  
+  for(Int_t i=0;i<fNCutSteps[firstVar];i++){
+    for(Int_t j=0;j<fNCutSteps[secondVar];j++){
+      index[firstVar]=i;
+      index[secondVar]=j;
+      Float_t cont=GetElement(index,ptbin)/norm;
+      Int_t bin1=i+1;
+      if(!fGreaterThan[firstVar]) bin1=fNCutSteps[firstVar]-i;
+      Int_t bin2=j+1;
+      if(!fGreaterThan[secondVar]) bin2=fNCutSteps[secondVar]-j;
+      h2->SetBinContent(bin1,bin2,cont);
+    }
+  }
+  return h2;
+}
+//_____________________________________________________________________________ 
+void  AliMultiDimVector::GetIntegrationLimits(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
+  // computes bin limits for integrating the AliMultiDimVector
+  minbin=0;
+  maxbin=0;
+  if(iVar<fNVariables){
+    minbin=iCell;
+    maxbin=fNCutSteps[iVar]-1;
+  }
+}
+//_____________________________________________________________________________ 
+void  AliMultiDimVector::GetFillRange(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const {
+  // computes range of cells passing the cuts for FillAndIntegrate
+  minbin=0;
+  maxbin=0;
+  if(iVar<fNVariables){
+    minbin=0; // bin 0 corresponds to loose cuts
+    maxbin=iCell;
+  }
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::Integrate(){
+  // integrates the matrix
+  if(fIsIntegrated){
+    AliError("MultiDimVector already integrated");
+    return;
+  }
+  TArrayF integral(fNTotCells);
+  for(ULong64_t i=0;i<fNTotCells;i++) integral[i]=CountsAboveCell(i);
+  for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]= integral[i];
+  fIsIntegrated=kTRUE;
+}
+//_____________________________________________________________________________ 
+Float_t AliMultiDimVector::CountsAboveCell(ULong64_t globadd) const{
+  // integrates the counts of cells above cell with address globadd
+  Int_t ind[fgkMaxNVariables];
+  Int_t ptbin;
+  GetIndicesFromGlobalAddress(globadd,ind,ptbin);
+  for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
+  Int_t mink[fgkMaxNVariables];
+  Int_t maxk[fgkMaxNVariables];
+  for(Int_t i=0;i<fgkMaxNVariables;i++){
+    GetIntegrationLimits(i,ind[i],mink[i],maxk[i]);
+  }
+  Float_t sumcont=0.;
+  for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
+    for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
+      for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
+       for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
+         for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
+           for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
+             for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
+               for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
+                 for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
+                   for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
+                     Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
+                     sumcont+=GetElement(currentBin,ptbin);
+                   }
+                 }
+               }
+             }
+           }
+         }
+       }
+      }
+    }
+  }
+  return sumcont;
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::Fill(Float_t* values, Int_t ptbin){
+  // fills the cells of AliMultiDimVector corresponding to values
+  if(fIsIntegrated){
+    AliError("MultiDimVector already integrated -- Use FillAndIntegrate");
+    return;
+  }
+  Int_t ind[fgkMaxNVariables];
+  Bool_t retcode=GetIndicesFromValues(values,ind);
+  for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
+  if(retcode) IncrementElement(ind,ptbin);
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::FillAndIntegrate(Float_t* values, Int_t ptbin){
+  // fills the cells of AliMultiDimVector passing the cuts
+  // The number of nested loops must match fgkMaxNVariables!!!!
+  fIsIntegrated=kTRUE;
+  Int_t ind[fgkMaxNVariables];
+  Bool_t retcode=GetIndicesFromValues(values,ind);
+  if(!retcode) return;
+  for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
+  Int_t mink[fgkMaxNVariables];
+  Int_t maxk[fgkMaxNVariables];
+  for(Int_t i=0;i<fgkMaxNVariables;i++){
+    GetFillRange(i,ind[i],mink[i],maxk[i]);
+  }
+  for(Int_t k0=mink[0]; k0<=maxk[0]; k0++){
+    for(Int_t k1=mink[1]; k1<=maxk[1]; k1++){
+      for(Int_t k2=mink[2]; k2<=maxk[2]; k2++){
+       for(Int_t k3=mink[3]; k3<=maxk[3]; k3++){
+         for(Int_t k4=mink[4]; k4<=maxk[4]; k4++){
+           for(Int_t k5=mink[5]; k5<=maxk[5]; k5++){
+             for(Int_t k6=mink[6]; k6<=maxk[6]; k6++){
+               for(Int_t k7=mink[7]; k7<=maxk[7]; k7++){
+                 for(Int_t k8=mink[8]; k8<=maxk[8]; k8++){
+                   for(Int_t k9=mink[9]; k9<=maxk[9]; k9++){
+                     Int_t currentBin[fgkMaxNVariables]={k0,k1,k2,k3,k4,k5,k6,k7,k8,k9};
+                     IncrementElement(currentBin,ptbin);
+                   }
+                 }
+               }
+             }
+           }
+         }
+       }
+      }
+    }
+  }
+
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::SuppressZeroBKGEffect(AliMultiDimVector* mvBKG){
+  // Sets to zero elements for which mvBKG=0
+  for(ULong64_t i=0;i<fNTotCells;i++)
+    if(mvBKG->GetElement(i)==0) fVett.AddAt(0,i);
+}
+//_____________________________________________________________________________ 
+AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBin){
+  // sums the elements of pt bins between firstBin and lastBin
+  if(firstBin<0 || lastBin>=fNPtBins || firstBin>=lastBin){
+    AliError("Bad numbers of Pt bins to be shrinked");
+    return 0;
+  }
+  Int_t nofcells[fgkMaxNVariables];
+  Float_t loosecuts[fgkMaxNVariables];
+  Float_t tightcuts[fgkMaxNVariables];
+  TString axisTitles[fgkMaxNVariables];
+  for(Int_t i=0;i<fNVariables;i++){
+    nofcells[i]=fNCutSteps[i];
+    if(fGreaterThan[i]){
+      loosecuts[i]=fMinLimits[i];
+      tightcuts[i]=fMaxLimits[i];
+    }else{
+      loosecuts[i]=fMaxLimits[i];
+      tightcuts[i]=fMinLimits[i];
+    }
+    axisTitles[i]=fAxisTitles[i];
+  }
+  Int_t newNptbins=fNPtBins-(lastBin-firstBin);
+  AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),fNVariables,newNptbins,nofcells,loosecuts,tightcuts,axisTitles);
+  
+  ULong64_t nOfPointsPerPtbin=fNTotCells/fNPtBins;
+  ULong64_t addressOld,addressNew;
+  Int_t npb,opb;
+  for(npb=0;npb<firstBin;npb++){
+    opb=npb;
+    for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
+      addressOld=opb+k*fNPtBins;
+      addressNew=npb+k*newNptbins;
+      shrinkedMV->SetElement(addressNew,fVett[addressOld]);
+    }
+  }
+  npb=firstBin;
+  for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
+    Float_t summedValue=0.;
+    for(opb=firstBin;opb<=lastBin;opb++){
+      addressOld=opb+k*fNPtBins;
+      summedValue+=fVett[addressOld];
+    }
+    addressNew=npb+k*newNptbins;
+    shrinkedMV->SetElement(addressNew,summedValue);
+  }
+  for(npb=firstBin+1;npb<newNptbins;npb++){
+    opb=npb+(lastBin-firstBin);
+    for(ULong64_t k=0;k<nOfPointsPerPtbin;k++){
+      addressOld=opb+k*fNPtBins;
+      addressNew=npb+k*newNptbins;
+      shrinkedMV->SetElement(addressNew,fVett[addressOld]);
+    }
+  }
+  return shrinkedMV;
+}
diff --git a/PWG3/vertexingHF/AliMultiDimVector.h b/PWG3/vertexingHF/AliMultiDimVector.h
new file mode 100644 (file)
index 0000000..ae4570f
--- /dev/null
@@ -0,0 +1,125 @@
+#ifndef ALIMULTIDIMVECTOR_H
+#define ALIMULTIDIMVECTOR_H
+
+/* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Class to store number of signal and background candidates     //
+// in bins of cut variables                                      //
+// Origin:       Elena Bruna (bruna@to.infn.it)                  //
+// Updated:      Sergey Senyukov (senyukov@to.infn.it)           //
+// Last updated: Francesco Prino (prino@to.infn.it)              //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+#include "TArrayF.h"
+#include "TArrayI.h"
+#include "TNamed.h"
+#include "TH2.h"
+#include "TString.h"
+
+class AliMultiDimVector :  public TNamed{
+
+ public:
+  AliMultiDimVector();
+  AliMultiDimVector(const AliMultiDimVector &mv);
+  AliMultiDimVector(const char *name, const char *title, const Int_t npars, Int_t nptbins, Int_t *nofcells, Float_t *loosecuts, Float_t *tightcuts, TString *axisTitles);
+  virtual ~AliMultiDimVector(){};
+
+  ULong64_t GetNTotCells()            const {return fNTotCells;}
+  Int_t     GetNVariables()           const {return fNVariables;}
+  Int_t     GetNPtBins()              const {return fNPtBins;}
+  Int_t     GetNCutSteps(Int_t iVar)  const {return fNCutSteps[iVar];}
+  Float_t   GetMinLimit(Int_t iVar)   const {return fMinLimits[iVar];}
+  Float_t   GetMaxLimit(Int_t iVar)   const {return fMaxLimits[iVar];}
+  Float_t   GetCutStep(Int_t iVar)    const {return (fMaxLimits[iVar]-fMinLimits[iVar])/(Float_t)fNCutSteps[iVar];}
+  TString   GetAxisTitle(Int_t iVar)  const {return fAxisTitles[iVar];}
+  Bool_t    IsIntegrated()            const {return fIsIntegrated;}
+
+  void CopyStructure(const AliMultiDimVector* mv);
+
+  Float_t   GetCutValue(Int_t iVar, Int_t iCell) const{
+    if(fGreaterThan[iVar]) return fMinLimits[iVar]+(Float_t)iCell*GetCutStep(iVar);
+    else return fMaxLimits[iVar]-(Float_t)iCell*GetCutStep(iVar);
+  }
+  Float_t   GetElement(ULong64_t globadd) const {return fVett[globadd];}
+  Float_t   GetElement(Int_t *ind, Int_t ptbin) const {
+    ULong64_t elem=GetGlobalAddressFromIndices(ind,ptbin);
+    return fVett[elem];
+  }
+  void      GetEntireMultiDimVector(Float_t *vett) const {
+    for(ULong64_t i=0; i<fNTotCells; i++) vett[i]=fVett[i];
+  }
+
+  Bool_t    GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const;
+  ULong64_t GetGlobalAddressFromIndices(Int_t *ind, Int_t ptbin) const;
+  Bool_t    GetIndicesFromValues(Float_t *values, Int_t *ind) const;
+  ULong64_t GetGlobalAddressFromValues(Float_t *values, Int_t ptbin) const;
+  Bool_t    GetCutValuesFromGlobalAddress(ULong64_t globadd, Float_t *cuts, Int_t &ptbin) const;
+
+  void SetElement(ULong64_t globadd,Float_t val) {fVett[globadd]=val;}
+  void SetElement(Int_t *ind, Int_t ptbin, Float_t val){
+    ULong64_t elem=GetGlobalAddressFromIndices(ind,ptbin);
+    fVett[elem]=val;
+  }
+  void IncrementElement(Int_t *ind, Int_t ptbin){
+    SetElement(ind,ptbin,GetElement(ind,ptbin)+1);
+  }
+  void IncrementElement(ULong64_t globadd){
+    SetElement(globadd,GetElement(globadd)+1.);
+  }
+
+  void Fill(Float_t* values, Int_t ptbin);
+  void FillAndIntegrate(Float_t* values, Int_t ptbin);
+  void Integrate();
+
+  void Reset(){
+    for(ULong64_t i=0; i<fNTotCells; i++) fVett[i]=0.;
+  }
+  void MultiplyBy(Float_t factor);
+  void Multiply(AliMultiDimVector* mv,Float_t factor);
+  void Multiply(AliMultiDimVector* mv1, AliMultiDimVector* mv2);
+  void Add(AliMultiDimVector* mv);
+  void Sum(AliMultiDimVector* mv1, AliMultiDimVector* mv2);
+  void LinearComb(AliMultiDimVector* mv1, Float_t norm1, AliMultiDimVector* mv2, Float_t norm2);
+  void DivideBy(AliMultiDimVector* mv);
+  void Divide(AliMultiDimVector* mv1, AliMultiDimVector* mv2);
+  void Sqrt();
+  void Sqrt(AliMultiDimVector* mv);
+  
+  void FindMaximum(Float_t& max_value, Int_t *ind, Int_t ptbin); 
+
+  TH2F*  Project(Int_t firstVar, Int_t secondVar, Int_t* fixedVars, Int_t ptbin, Float_t norm=1.);
+
+  void SuppressZeroBKGEffect(AliMultiDimVector* BKG);
+   AliMultiDimVector* ShrinkPtBins(Int_t firstBin, Int_t lastBin);
+
+ protected:
+  void GetIntegrationLimits(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const;
+  void GetFillRange(Int_t iVar, Int_t iCell, Int_t& minbin, Int_t& maxbin) const;
+  Bool_t    GetGreaterThan(Int_t iVar) const {return fGreaterThan[iVar];}
+  Float_t   CountsAboveCell(ULong64_t globadd) const;
+
+ private:
+  static const Int_t fgkMaxNVariables=10;  // max. n. of selection variables
+
+  Int_t     fNVariables;                   // n. of selection variables
+  Int_t     fNPtBins;                      // n. of pt bins
+  Int_t     fNCutSteps[fgkMaxNVariables];  // n. of cut step for each variable
+  Float_t   fMinLimits[fgkMaxNVariables];  // lower cut value for each variable
+  Float_t   fMaxLimits[fgkMaxNVariables];  // higher cut value for each variable
+  Bool_t    fGreaterThan[fgkMaxNVariables];// sign of the cut (> or <)
+  TString   fAxisTitles[fgkMaxNVariables]; // titles for variables
+  TArrayF   fVett;                   // array with n. of candidates vs. cuts
+  ULong64_t fNTotCells;              // total number of matrix elements
+  Bool_t    fIsIntegrated;           // flag for integrated matrix 
+
+  ClassDef(AliMultiDimVector,1); // a multi-dimensional vector class
+
+};
+
+#endif
diff --git a/PWG3/vertexingHF/AliSignificanceCalculator.cxx b/PWG3/vertexingHF/AliSignificanceCalculator.cxx
new file mode 100644 (file)
index 0000000..628a701
--- /dev/null
@@ -0,0 +1,219 @@
+/**************************************************************************
+ * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Implementation of the class to calculate statistical          //
+// significance from AliMultiVeector objects with signal and     //
+// background counts vs. cut values                              //
+// Origin: Francesco Prino (prino@to.infn.it)                    //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+#include "AliSignificanceCalculator.h"
+#include "TMath.h"
+#include "AliLog.h"
+
+ClassImp(AliSignificanceCalculator)
+//___________________________________________________________________________
+AliSignificanceCalculator::AliSignificanceCalculator():
+fSignal(0),
+fErrSquareSignal(0),
+fBackground(0),
+fErrSquareBackground(0),
+fSignificance(0),
+fErrSignificance(0),
+fNormSig(1.),  
+fNormBkg(1.)
+{
+  // default constructor
+  fSignal=new AliMultiDimVector();
+  fBackground=new AliMultiDimVector();
+}
+//___________________________________________________________________________
+AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, Float_t normsig, Float_t normbkg):
+fSignal(sig),
+fErrSquareSignal(0),
+fBackground(bkg),
+fErrSquareBackground(0),
+fSignificance(0),
+fErrSignificance(0),
+fNormSig(normsig),
+fNormBkg(normbkg)
+{
+  // standard constructor
+  if(fSignal && fBackground) CalculateSignificance();
+}
+//___________________________________________________________________________
+AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, AliMultiDimVector* err2sig, AliMultiDimVector* err2bkg, Float_t normsig, Float_t normbkg):
+fSignal(sig),
+fErrSquareSignal(err2sig),
+fBackground(bkg),
+fErrSquareBackground(err2bkg),
+fSignificance(0),
+fErrSignificance(0),
+fNormSig(normsig),
+fNormBkg(normbkg)
+{
+  // standard constructor
+  if(fSignal && fBackground) CalculateSignificance();
+}
+//___________________________________________________________________________
+AliSignificanceCalculator::~AliSignificanceCalculator(){
+  // destructor
+  if(fSignal) delete fSignal;
+  if(fBackground) delete fBackground;
+  if(fSignificance) delete fSignificance;
+  if(fErrSignificance) delete fErrSignificance;
+}
+//___________________________________________________________________________
+Bool_t AliSignificanceCalculator::Check() const {
+  // checks AliMultiDimVector dimension and normalization
+  if(fSignal==0 || fBackground==0) return kFALSE;
+  if(fNormSig==0. || fNormBkg==0.) return kFALSE;
+  if(fSignal->GetNTotCells() != fBackground->GetNTotCells()) return kFALSE;
+  return kTRUE;
+}
+//___________________________________________________________________________
+void AliSignificanceCalculator::CalculateSignificance(){
+  // calculates significance and its error
+  if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
+  
+  if(fSignificance) delete fSignificance;
+  if(fErrSignificance) delete fErrSignificance;
+  fSignificance=new AliMultiDimVector();
+  fSignificance->CopyStructure(fSignal);
+  fErrSignificance=new AliMultiDimVector();
+  fErrSignificance->CopyStructure(fSignal);
+
+  for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
+    if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
+      Float_t s=fSignal->GetElement(i)*fNormSig;
+      Float_t b=fBackground->GetElement(i)*fNormBkg;
+      Float_t signif=0.;
+      Float_t errsig=0.;
+      if((s+b)>0){ 
+       signif=s/TMath::Sqrt(s+b);
+       Float_t errs,errb;
+       if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
+       else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
+       if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
+       else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
+       Float_t dsigds=(s+2*b)/2./(s+b)/TMath::Sqrt(s+b);
+       Float_t dsigdb=-s/2./(s+b)/TMath::Sqrt(s+b);
+       errsig=TMath::Sqrt(dsigds*dsigds*errs*errs+dsigdb*dsigdb*errb*errb);
+      }
+      fSignificance->SetElement(i,signif);
+      fErrSignificance->SetElement(i,errsig);
+    }
+  }
+  fSignificance->SetNameTitle("Significance","Significance");
+  fErrSignificance->SetNameTitle("ErrorOnSignificance","ErrorOnSignificance");
+}
+//___________________________________________________________________________
+AliMultiDimVector* AliSignificanceCalculator::CalculatePurity() const {
+  // calculates purity 
+  if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
+  
+  AliMultiDimVector* purity=new AliMultiDimVector();
+  purity->CopyStructure(fSignal);
+  for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
+    if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
+      Float_t s=fSignal->GetElement(i)*fNormSig;
+      Float_t b=fBackground->GetElement(i)*fNormBkg;
+      Float_t pur=0.;
+      if((s+b)>0) pur=s/(s+b);
+      purity->SetElement(i,pur);
+    }
+  }
+  purity->SetNameTitle("Purity","Purity");
+  return purity;
+}
+//___________________________________________________________________________
+AliMultiDimVector* AliSignificanceCalculator::CalculatePurityError() const {
+  // calculates error on purity 
+  if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
+  
+  AliMultiDimVector* epurity=new AliMultiDimVector();
+  epurity->CopyStructure(fSignal);
+  for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
+    if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
+      Float_t s=fSignal->GetElement(i)*fNormSig;
+      Float_t b=fBackground->GetElement(i)*fNormBkg;
+      Float_t epur=0.;
+      if((s+b)>0){
+       Float_t errs,errb;
+       if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
+       else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
+       if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
+       else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
+       Float_t dpurds=b/(s+b)/(s+b);
+       Float_t dpurdb=-s/(s+b)/(s+b);
+       epur=TMath::Sqrt(dpurds*dpurds*errs*errs+dpurdb*dpurdb*errb*errb);
+      }
+      epurity->SetElement(i,epur);
+    }
+  }
+  epurity->SetNameTitle("ErrorOnPurity","ErrorOnPurity");
+  return epurity;
+}
+//___________________________________________________________________________
+AliMultiDimVector* AliSignificanceCalculator::CalculateSOverB() const {
+  // Signal over Background
+  if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
+  
+  AliMultiDimVector* sob=new AliMultiDimVector();
+  sob->CopyStructure(fSignal);
+  for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
+    if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
+      Float_t s=fSignal->GetElement(i)*fNormSig;
+      Float_t b=fBackground->GetElement(i)*fNormBkg;
+      Float_t soverb=0.;
+      if(b>0) soverb=s/b;
+      sob->SetElement(i,soverb);
+    }
+  }
+  sob->SetNameTitle("SoverB","SoverB");
+  return sob;
+}
+//___________________________________________________________________________
+AliMultiDimVector* AliSignificanceCalculator::CalculateSOverBError() const {
+  // Error on Signal over Background
+  if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
+  
+  AliMultiDimVector* esob=new AliMultiDimVector();
+  esob->CopyStructure(fSignal);
+  for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
+    if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
+      Float_t s=fSignal->GetElement(i)*fNormSig;
+      Float_t b=fBackground->GetElement(i)*fNormBkg;
+      Float_t esoverb=0.;
+      if(b>0){
+       Float_t soverb=s/b;
+       Float_t errs,errb;
+       if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
+       else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
+       if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
+       else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
+       esoverb=soverb*TMath::Sqrt(errs*errs/s/s+errb*errb/b/b);
+      }
+      esob->SetElement(i,esoverb);
+    }
+  }
+  esob->SetNameTitle("ErrorOnSoverB","ErrorOnSoverB");
+  return esob;
+}
diff --git a/PWG3/vertexingHF/AliSignificanceCalculator.h b/PWG3/vertexingHF/AliSignificanceCalculator.h
new file mode 100644 (file)
index 0000000..82e6b58
--- /dev/null
@@ -0,0 +1,93 @@
+#ifndef ALISIGNIFICANCECALCULATOR_H
+#define ALISIGNIFICANCECALCULATOR_H
+
+/* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Class to calculate the statistical significance from          //
+// AliMultiVeector objects for signal and background             //
+// Origin: Francesco Prino (prino@to.infn.it)                    //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+#include "AliMultiDimVector.h"
+#include "TObject.h"
+
+class AliSignificanceCalculator : public TObject{
+ public:
+  AliSignificanceCalculator();
+  AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, Float_t normsig=1., Float_t normbkg=1.);
+  AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, AliMultiDimVector* err2sig, AliMultiDimVector* err2bkg, Float_t normsig=1., Float_t normbkg=1.);
+
+  ~AliSignificanceCalculator();
+
+  void SetSignal(AliMultiDimVector* sig, Float_t norm=1.){
+    if(fSignal) delete fSignal;
+    fSignal=sig;
+    fNormSig=norm;
+    if(fSignal && fBackground) CalculateSignificance();
+  }
+  void SetBackground(AliMultiDimVector* bac, Float_t norm=1.){
+    if(fBackground) delete fBackground;
+    fBackground=bac;
+    fNormBkg=norm;
+    if(fSignal && fBackground) CalculateSignificance();
+  }
+  void SetErrSquareSignal(AliMultiDimVector* err2sig, Float_t norm=1.){
+    if(fErrSquareSignal) delete fErrSquareSignal;
+    fErrSquareSignal=err2sig;
+    fNormSig=norm;
+    if(fSignal && fBackground) CalculateSignificance();
+  }
+  void SetErrSquareBackground(AliMultiDimVector* err2bkg, Float_t norm=1.){
+    if(fErrSquareBackground) delete fErrSquareBackground;
+    fErrSquareBackground=err2bkg;
+    fNormBkg=norm;
+    if(fSignal && fBackground) CalculateSignificance();
+  }
+  
+  void SetNormalizations(Float_t normSig, Float_t normBkg){
+    fNormSig=normSig;
+    fNormBkg=normBkg;
+    if(fSignal && fBackground) CalculateSignificance();
+  }
+
+  AliMultiDimVector* GetSignal() const {return fSignal;}
+  AliMultiDimVector* GetBackground() const {return fBackground;}
+  AliMultiDimVector* GetSignificance() const {return fSignificance;}
+  AliMultiDimVector* GetSignificanceError() const {return fErrSignificance;}
+
+  void CalculateSignificance();
+  Float_t GetMaxSignificance(Int_t* cutIndices, Int_t ptbin) const{
+    Float_t sigMax=0;
+    if(fSignificance) fSignificance->FindMaximum(sigMax,cutIndices,ptbin);
+    return sigMax;
+  }
+  AliMultiDimVector* CalculatePurity() const;
+  AliMultiDimVector* CalculatePurityError() const;
+  AliMultiDimVector* CalculateSOverB() const;
+  AliMultiDimVector* CalculateSOverBError() const;
+
+ private:
+  Bool_t Check() const;
+  AliSignificanceCalculator(const AliSignificanceCalculator& c);
+  AliSignificanceCalculator& operator=(const AliSignificanceCalculator& c);
+
+  AliMultiDimVector* fSignal;              // signal matrix
+  AliMultiDimVector* fErrSquareSignal;     // matrix with err^2 for signal
+  AliMultiDimVector* fBackground;          // background matrix
+  AliMultiDimVector* fErrSquareBackground; // matrix with err^2 for background
+  AliMultiDimVector* fSignificance;        // significance matrix
+  AliMultiDimVector* fErrSignificance;     // matrix with error on significance
+  Float_t fNormSig;                        // signal normalization
+  Float_t fNormBkg;                        // background normalization
+
+  ClassDef(AliSignificanceCalculator,1); // class to compute and maximise significance
+
+};
+
+#endif
diff --git a/PWG3/vertexingHF/TestMultiVector.C b/PWG3/vertexingHF/TestMultiVector.C
new file mode 100644 (file)
index 0000000..c6910f2
--- /dev/null
@@ -0,0 +1,132 @@
+void TestMultiVector(){
+
+  // Example of usage of AliMultiDimVector and AliSignificanceCalculator classes
+
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libAOD");
+  gSystem->Load("libCORRFW");
+  gSystem->Load("libPWG3");
+  gSystem->Load("libPWG3vertexingHF");
+
+  Int_t nptbins=4;
+  const Int_t npars=3;
+  Int_t nofcells[npars]={20,20,20};
+  Float_t looses[npars]={0.,700.,0.8};
+  Float_t tights[npars]={1500.,0.,1.};
+  TString AxisTitle[npars];
+  AxisTitle[0]="DistPS (micron)";
+  AxisTitle[1]="TrackDispersion";
+  AxisTitle[2]="CosPoint";
+
+  AliMultiDimVector* mvsig=new AliMultiDimVector("Signal","Signal",npars,nptbins,nofcells,looses,tights,AxisTitle);
+  AliMultiDimVector* mvsig2=new AliMultiDimVector("Signal","Signal",npars,nptbins,nofcells,looses,tights,AxisTitle);
+  AliMultiDimVector* mvbkg=new AliMultiDimVector("Background","Background",npars,nptbins,nofcells,looses,tights,AxisTitle);
+
+
+  TF1* dsig=new TF1("dsig","exp(-x/310.)",0.,5000.);
+  TF1* dbak=new TF1("dbak","exp(-x/50)",0.,5000.);
+  gStyle->SetOptStat(2);
+  for(Int_t isignev=0;isignev<5000;isignev++){
+    Int_t ptbin=gRandom->Integer(nptbins);
+    Float_t dists=dsig->GetRandom();
+    Float_t distb=dbak->GetRandom();
+    Float_t cpas=1-TMath::Abs(gRandom->Gaus(0.,0.02));
+    Float_t cpab=0.8+gRandom->Rndm()*0.2;
+    Float_t sigverts=gRandom->Gaus(200.,25.);
+    Float_t sigvertb=gRandom->Gaus(250.,25.);
+    Float_t vals[npars]={dists,sigverts,cpas};
+    Float_t valb[npars]={distb,sigvertb,cpab};
+    mvsig->FillAndIntegrate(vals,ptbin);
+    mvsig2->Fill(vals,ptbin);  // alternative way of filling
+    mvbkg->FillAndIntegrate(valb,ptbin);
+  }
+  mvsig2->Integrate(); // mvsig2 is now equal to mvsig
+  
+  // Merge Pt bins
+  AliMultiDimVector* mvsigallpt=mvsig->ShrinkPtBins(0,3);
+  AliMultiDimVector* mvbkgallpt=mvbkg->ShrinkPtBins(0,3);
+  
+  //caluclate significance
+  AliSignificanceCalculator* cal=new AliSignificanceCalculator(mvsigallpt,mvbkgallpt,1.,5.);
+  AliMultiDimVector* mvsts=cal->GetSignificance();
+  AliMultiDimVector* mvess=cal->GetSignificanceError();
+  AliMultiDimVector* mvsob=cal->CalculateSOverB();
+  AliMultiDimVector* mvesob=cal->CalculateSOverBError();
+  AliMultiDimVector* mvpur=cal->CalculatePurity();
+  AliMultiDimVector* mvepur=cal->CalculatePurityError();
+  Int_t fixed[3]={0,0,0};
+  gStyle->SetPalette(1);
+
+  Int_t maxInd[3];
+  Float_t sigMax=cal->GetMaxSignificance(maxInd,0);
+  Float_t cut0=mvsigallpt->GetCutValue(0,maxInd[0]);
+  Float_t cut1=mvsigallpt->GetCutValue(1,maxInd[1]);
+  Float_t cut2=mvsigallpt->GetCutValue(2,maxInd[2]);
+
+  printf("=========== Pt Integrated ==============\n");
+  printf("Maximum of significance found in bin %d %d %d\n",maxInd[0],maxInd[1],maxInd[2]);
+  printf("                                 cuts= %f %f %f\n",cut0,cut1,cut2);
+  printf("Significance = %f +- %f\n",mvsts->GetElement(maxInd,0),mvess->GetElement(maxInd,0));
+  printf("S/B          = %f +- %f\n",mvsob->GetElement(maxInd,0),mvesob->GetElement(maxInd,0));
+  printf("Purity       = %f +- %f\n",mvpur->GetElement(maxInd,0),mvepur->GetElement(maxInd,0));
+
+
+  // 2D projections
+  TH2F* hsig1 = mvsigallpt->Project(0,2,fixed,0);
+  TH2F* hbkg1 = mvbkgallpt->Project(0,2,fixed,0);
+  TH2F* hsts1 = mvsts->Project(0,2,fixed,0);
+  TH2F* hess1 = mvess->Project(0,2,fixed,0);
+  TH2F* hpur1 = mvpur->Project(0,2,fixed,0);
+  TH2F* hepur1 = mvepur->Project(0,2,fixed,0);
+  TH2F* hsob1 = mvsob->Project(0,2,fixed,0);
+  TH2F* hesob1 = mvesob->Project(0,2,fixed,0);
+
+  TH2F* hsig2 = mvsigallpt->Project(1,2,fixed,0);
+  TH2F* hbkg2 = mvbkgallpt->Project(1,2,fixed,0);
+  TH2F* hsts2 = mvsts->Project(1,2,fixed,0);
+  TH2F* hess2 = mvess->Project(1,2,fixed,0);
+  TH2F* hpur2 = mvpur->Project(1,2,fixed,0);
+  TH2F* hepur2 = mvepur->Project(1,2,fixed,0);
+  TH2F* hsob2 = mvsob->Project(1,2,fixed,0);
+  TH2F* hesob2 = mvesob->Project(1,2,fixed,0);
+
+  TCanvas* c1=new TCanvas("c1","Var 0 vs. 2",1000,800);
+  c1->Divide(4,2);
+  c1->cd(1);
+  hsig1->Draw("colz");
+  c1->cd(2);
+  hbkg1->Draw("colz");
+  c1->cd(3);
+  hsts1->Draw("colz");
+  c1->cd(4);
+  hess1->Draw("colz");
+  c1->cd(5);
+  hsob1->Draw("colz");
+  c1->cd(6);
+  hesob1->Draw("colz");
+  c1->cd(7);
+  hpur1->Draw("colz");
+  c1->cd(8);
+  hepur1->Draw("colz");
+
+  TCanvas* c2=new TCanvas("c2","Var 1 vs. 2",1000,800);
+  c2->Divide(4,2);
+  c2->cd(1);
+  hsig2->Draw("colz");
+  c2->cd(2);
+  hbkg2->Draw("colz");
+  c2->cd(3);
+  hsts2->Draw("colz");
+  c2->cd(4);
+  hess2->Draw("colz");
+  c2->cd(5);
+  hsob2->Draw("colz");
+  c2->cd(6);
+  hesob2->Draw("colz");
+  c2->cd(7);
+  hpur2->Draw("colz");
+  c2->cd(8);
+  hepur2->Draw("colz");
+
+}