]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliMultiDimVector.cxx
Update (Zaida)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliMultiDimVector.cxx
index b6d61be32b1f260bcf5aea63f90631b6d3ea0644..4d32ea33830fcad13209484cb0748c0e3d137669 100644 (file)
@@ -21,7 +21,8 @@
 // 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)              //
+//               Francesco Prino (prino@to.infn.it)              //
+// Last Updated: Giacomo Ortona (ortona@to.infn.it)              //
 //                                                               //
 ///////////////////////////////////////////////////////////////////
 
@@ -29,7 +30,6 @@
 #include "TH2.h"
 #include "AliMultiDimVector.h"
 #include "AliLog.h"
-#include "TMath.h"
 #include "TString.h"
 
 ClassImp(AliMultiDimVector)
@@ -44,7 +44,7 @@ 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),
+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),
 fNVariables(npars),
 fNPtBins(nptbins),
 fVett(0),
@@ -55,7 +55,18 @@ fIsIntegrated(0){
   for(Int_t i=0;i<fNVariables;i++){
     ntot*=nofcells[i];
     fNCutSteps[i]=nofcells[i];
-    if(loosecuts[i] <= tightcuts[i]){
+    if(loosecuts[i] == tightcuts[i]){
+      if(loosecuts[i]!=0){
+       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]);
+       fMinLimits[i]=tightcuts[i]-0.1*tightcuts[i];
+       fMaxLimits[i]=tightcuts[i];
+      }else{
+       fMinLimits[i]=0;
+       fMaxLimits[i]=0.0001;
+      }
+       fGreaterThan[i]=kTRUE;
+    }
+    if(loosecuts[i] < tightcuts[i]){
       fMinLimits[i]=loosecuts[i];
       fMaxLimits[i]=tightcuts[i];
       fGreaterThan[i]=kTRUE;
@@ -66,8 +77,10 @@ fIsIntegrated(0){
     }
     fAxisTitles[i]=axisTitles[i].Data();
   }
-  fNTotCells=ntot*fNPtBins;
+  fNTotCells=ntot*fNPtBins;  
   fVett.Set(fNTotCells); 
+  for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=ptlimits[ipt];
+  for(Int_t ipt=fNPtBins+1;ipt<fgkMaxNPtBins+1;ipt++) fPtLimits[ipt]=999.;
   for (ULong64_t j=0;j<fNTotCells;j++) fVett.AddAt(0,j);
 }
 //___________________________________________________________________________
@@ -87,6 +100,8 @@ fIsIntegrated(mv.fIsIntegrated)
     fAxisTitles[i]=mv.GetAxisTitle(i);
   }
   fVett.Set(fNTotCells); 
+  
+  for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv.GetPtLimit(ipt);
   for(ULong64_t i=0;i<fNTotCells;i++) fVett[i]=mv.GetElement(i);
 }
 //___________________________________________________________________________
@@ -103,8 +118,8 @@ void AliMultiDimVector::CopyStructure(const AliMultiDimVector* mv){
     fGreaterThan[i]=mv->GetGreaterThan(i);
     fAxisTitles[i]=mv->GetAxisTitle(i);
   }
-  fVett.Set(fNTotCells);
-  
+  for(Int_t ipt=0;ipt<fNPtBins+1;ipt++) fPtLimits[ipt]=mv->GetPtLimit(ipt);
+  fVett.Set(fNTotCells);  
 }
 //______________________________________________________________________
 Bool_t AliMultiDimVector::GetIndicesFromGlobalAddress(ULong64_t globadd, Int_t *ind, Int_t &ptbin) const {
@@ -134,7 +149,7 @@ Bool_t AliMultiDimVector::GetCutValuesFromGlobalAddress(ULong64_t globadd, Float
   return kTRUE;
 }
 //______________________________________________________________________
-ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(Int_t *ind, Int_t ptbin) const {
+ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(const 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;
@@ -159,7 +174,7 @@ ULong64_t AliMultiDimVector::GetGlobalAddressFromIndices(Int_t *ind, Int_t ptbin
   return elem;
 }
 //______________________________________________________________________
-Bool_t AliMultiDimVector::GetIndicesFromValues(Float_t *values, Int_t *ind) const {
+Bool_t AliMultiDimVector::GetIndicesFromValues(const 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]){ 
@@ -175,7 +190,7 @@ Bool_t AliMultiDimVector::GetIndicesFromValues(Float_t *values, Int_t *ind) cons
   return kTRUE;
 }
 //______________________________________________________________________
-ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(Float_t *values, Int_t ptbin) const {
+ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(const 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);
@@ -189,93 +204,93 @@ ULong64_t AliMultiDimVector::GetGlobalAddressFromValues(Float_t *values, Int_t p
 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)
+    if(fVett.At(i)>0.)
       fVett.AddAt(fVett.At(i)*factor,i);
     else fVett.AddAt(-1,i);
   }
   
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Multiply(AliMultiDimVector* mv,Float_t factor){
+void AliMultiDimVector::Multiply(const AliMultiDimVector* mv,Float_t factor){
   //  Sets AliMultiDimVector=mv*constant factor
   for(ULong64_t i=0;i<fNTotCells;i++){
-    if(mv->GetElement(i)!=-1)
+    if(mv->GetElement(i)>0.)
       fVett.AddAt(mv->GetElement(i)*factor,i);
     else fVett.AddAt(-1,i); 
   }
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Multiply(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+void AliMultiDimVector::Multiply(const AliMultiDimVector* mv1, const AliMultiDimVector* mv2){
   //  Sets AliMultiDimVector=mv1*mv2
   for(ULong64_t i=0;i<fNTotCells;i++){
-    if(mv1->GetElement(i)!=-1 && mv2->GetElement(i)!=-1)
+    if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
       fVett.AddAt(mv1->GetElement(i)*mv2->GetElement(i),i);
     else fVett.AddAt(-1,i); 
   }
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Add(AliMultiDimVector* mv){
+void AliMultiDimVector::Add(const 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)
+      if(mv->GetElement(i)>0. && fVett.At(i)>0.)
        fVett.AddAt(fVett.At(i)+mv->GetElement(i),i);
       else fVett.AddAt(-1,i); 
   }
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Sum(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+void AliMultiDimVector::Sum(const AliMultiDimVector* mv1, const 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)
+      if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
        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){
+void AliMultiDimVector::LinearComb(const AliMultiDimVector* mv1, Float_t norm1, const 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)
+      if(mv1->GetElement(i)>0. && mv2->GetElement(i)>0.)
        fVett.AddAt(norm1*mv1->GetElement(i)+norm2*mv2->GetElement(i),i); 
       else fVett.AddAt(-1,i); 
     }
   }
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::DivideBy(AliMultiDimVector* mv){
+void AliMultiDimVector::DivideBy(const 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)
+      if(mv->GetElement(i)!=0 &&mv->GetElement(i)>0. && fVett.At(i)>0.)
        fVett.AddAt(fVett.At(i)/mv->GetElement(i),i);
       else fVett.AddAt(-1,i);
   }
 
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Divide(AliMultiDimVector* mv1,AliMultiDimVector* mv2){
+void AliMultiDimVector::Divide(const AliMultiDimVector* mv1, const 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)
+      if(mv2->GetElement(i)!=0&& mv2->GetElement(i)>0.&& mv1->GetElement(i)>0.)
        {
          fVett.AddAt(mv1->GetElement(i)/mv2->GetElement(i),i);
        }
@@ -293,7 +308,7 @@ void AliMultiDimVector::Sqrt(){
   }
 }
 //_____________________________________________________________________________
-void AliMultiDimVector::Sqrt(AliMultiDimVector* mv){
+void AliMultiDimVector::Sqrt(const 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);
@@ -318,7 +333,7 @@ void AliMultiDimVector::FindMaximum(Float_t& maxValue, Int_t *ind , Int_t ptbin)
 }
 
 //_____________________________________________________________________________
-TH2F*  AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, Int_t* fixedVars, Int_t ptbin, Float_t norm){
+TH2F*  AliMultiDimVector::Project(Int_t firstVar, Int_t secondVar, const 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);
@@ -375,6 +390,49 @@ void AliMultiDimVector::Integrate(){
   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;
+}//_____________________________________________________________________________ 
+ULong64_t* AliMultiDimVector::GetGlobalAddressesAboveCuts(const Float_t *values, Int_t ptbin, Int_t& nVals) const{
+  // fills an array with global addresses of cells passing the cuts
+
+  Int_t ind[fgkMaxNVariables];
+  Bool_t retcode=GetIndicesFromValues(values,ind);
+  if(!retcode){ 
+    nVals=0;
+    return 0x0;
+  }
+  for(Int_t i=fNVariables; i<fgkMaxNVariables; i++) ind[i]=0;
+  Int_t mink[fgkMaxNVariables];
+  Int_t maxk[fgkMaxNVariables];
+  Int_t size=1;
+  for(Int_t i=0;i<fgkMaxNVariables;i++){
+    GetFillRange(i,ind[i],mink[i],maxk[i]);
+    size*=(maxk[i]-mink[i]+1);
+  }
+  ULong64_t* indexes=new ULong64_t[size];
+  nVals=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};
+                     indexes[nVals++]=GetGlobalAddressFromIndices(currentBin,ptbin);
+                   }
+                 }
+               }
+             }
+           }
+         }
+       }
+      }
+    }
+  }
+  return indexes;
 }
 //_____________________________________________________________________________ 
 Float_t AliMultiDimVector::CountsAboveCell(ULong64_t globadd) const{
@@ -464,10 +522,10 @@ void AliMultiDimVector::FillAndIntegrate(Float_t* values, Int_t ptbin){
 
 }
 //_____________________________________________________________________________ 
-void AliMultiDimVector::SuppressZeroBKGEffect(AliMultiDimVector* mvBKG){
+void AliMultiDimVector::SuppressZeroBKGEffect(const 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);
+    if(mvBKG->GetElement(i)<0.00000001) fVett.AddAt(0,i);
 }
 //_____________________________________________________________________________ 
 AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBin){
@@ -492,7 +550,10 @@ AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBi
     axisTitles[i]=fAxisTitles[i];
   }
   Int_t newNptbins=fNPtBins-(lastBin-firstBin);
-  AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),fNVariables,newNptbins,nofcells,loosecuts,tightcuts,axisTitles);
+  Float_t ptlimits[fgkMaxNPtBins+1];
+  for(Int_t ipt=0; ipt<=firstBin;ipt++) ptlimits[ipt]=fPtLimits[ipt];
+  for(Int_t ipt=firstBin+1; ipt<newNptbins+1;ipt++) ptlimits[ipt]=fPtLimits[ipt+(lastBin-firstBin)];
+  AliMultiDimVector* shrinkedMV=new AliMultiDimVector(GetName(),GetTitle(),newNptbins,ptlimits,fNVariables,nofcells,loosecuts,tightcuts,axisTitles);
   
   ULong64_t nOfPointsPerPtbin=fNTotCells/fNPtBins;
   ULong64_t addressOld,addressNew;
@@ -525,3 +586,39 @@ AliMultiDimVector* AliMultiDimVector:: ShrinkPtBins(Int_t firstBin, Int_t lastBi
   }
   return shrinkedMV;
 }
+//_____________________________________________________________________________ 
+void AliMultiDimVector::SetNewLimits(Float_t* loose,Float_t* tight){
+  for(Int_t i=0;i<fNVariables;i++){
+    if(loose[i] < tight[i]){
+      fMinLimits[i]=loose[i];
+      fMaxLimits[i]=tight[i];
+      fGreaterThan[i]=kTRUE;
+    }else{
+      fMinLimits[i]=tight[i];
+      fMaxLimits[i]=loose[i];
+      fGreaterThan[i]=kFALSE;
+    }
+  }
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::SwapLimits(Int_t ivar){
+  Float_t oldmin = fMinLimits[ivar];
+  fMinLimits[ivar] = fMaxLimits[ivar];
+  fMaxLimits[ivar] = oldmin;
+  if(fGreaterThan[ivar])fGreaterThan[ivar]=kFALSE;
+  else fGreaterThan[ivar]=kTRUE;
+}
+//_____________________________________________________________________________ 
+void AliMultiDimVector::PrintStatus(){
+  //
+  printf("Number of Pt bins       = %d\n",fNPtBins);
+  printf("Limits of Pt bins       = ");
+  for(Int_t ib=0;ib<fNPtBins+1;ib++) printf("%6.2f ",fPtLimits[ib]);
+  printf("\n");
+  printf("Number of cut variables = %d\n",fNVariables);
+  for(Int_t iv=0;iv<fNVariables;iv++){
+    printf("- Variable %d: %s\n",iv,fAxisTitles[iv].Data());
+    printf("    Nsteps= %d Rage = %6.2f %6.2f\n",
+          fNCutSteps[iv],fMinLimits[iv],fMaxLimits[iv]);
+  }
+}