--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+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");
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff