PWG3/vertexingHF/AliSignificanceCalculator.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: $ */
  17
  18 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // Implementation of the class to calculate statistical          //
  21 // significance from AliMultiVeector objects with signal and     //
  22 // background counts vs. cut values                              //
  23 // Origin: Francesco Prino (prino@to.infn.it)                    //
  24 //                                                               //
  25 ///////////////////////////////////////////////////////////////////
  26
  27 #include "AliSignificanceCalculator.h"
  28 #include "TMath.h"
  29 #include "AliLog.h"
  30
  31 ClassImp(AliSignificanceCalculator)
  32 //___________________________________________________________________________
  33 AliSignificanceCalculator::AliSignificanceCalculator():
  34 fSignal(0),
  35 fErrSquareSignal(0),
  36 fBackground(0),
  37 fErrSquareBackground(0),
  38 fSignificance(0),
  39 fErrSignificance(0),
  40 fNormSig(1.),
  41 fNormBkg(1.)
  42 {
  43   // default constructor
  44   fSignal=new AliMultiDimVector();
  45   fBackground=new AliMultiDimVector();
  46 }
  47 //___________________________________________________________________________
  48 AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, Float_t normsig, Float_t normbkg):
  49 fSignal(sig),
  50 fErrSquareSignal(0),
  51 fBackground(bkg),
  52 fErrSquareBackground(0),
  53 fSignificance(0),
  54 fErrSignificance(0),
  55 fNormSig(normsig),
  56 fNormBkg(normbkg)
  57 {
  58   // standard constructor
  59   if(fSignal && fBackground) CalculateSignificance();
  60 }
  61 //___________________________________________________________________________
  62 AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, AliMultiDimVector* err2sig, AliMultiDimVector* err2bkg, Float_t normsig, Float_t normbkg):
  63 fSignal(sig),
  64 fErrSquareSignal(err2sig),
  65 fBackground(bkg),
  66 fErrSquareBackground(err2bkg),
  67 fSignificance(0),
  68 fErrSignificance(0),
  69 fNormSig(normsig),
  70 fNormBkg(normbkg)
  71 {
  72   // standard constructor
  73   if(fSignal && fBackground) CalculateSignificance();
  74 }
  75 //___________________________________________________________________________
  76 AliSignificanceCalculator::~AliSignificanceCalculator(){
  77   // destructor
  78   if(fSignal) delete fSignal;
  79   if(fBackground) delete fBackground;
  80   if(fSignificance) delete fSignificance;
  81   if(fErrSignificance) delete fErrSignificance;
  82 }
  83 //___________________________________________________________________________
  84 Bool_t AliSignificanceCalculator::Check() const {
  85   // checks AliMultiDimVector dimension and normalization
  86   if(fSignal==0 || fBackground==0) return kFALSE;
  87   if(fNormSig==0. || fNormBkg==0.) return kFALSE;
  88   if(fSignal->GetNTotCells() != fBackground->GetNTotCells()) return kFALSE;
  89   return kTRUE;
  90 }
  91 //___________________________________________________________________________
  92 void AliSignificanceCalculator::CalculateSignificance(){
  93   // calculates significance and its error
  94   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
  95
  96   if(fSignificance) delete fSignificance;
  97   if(fErrSignificance) delete fErrSignificance;
  98   fSignificance=new AliMultiDimVector();
  99   fSignificance->CopyStructure(fSignal);
 100   fErrSignificance=new AliMultiDimVector();
 101   fErrSignificance->CopyStructure(fSignal);
 102
 103   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 104     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 105       Float_t s=fSignal->GetElement(i)*fNormSig;
 106       Float_t b=fBackground->GetElement(i)*fNormBkg;
 107       Float_t signif=0.;
 108       Float_t errsig=0.;
 109       if((s+b)>0){
 110         signif=s/TMath::Sqrt(s+b);
 111         Float_t errs,errb;
 112         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 113         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 114         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 115         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 116         Float_t dsigds=(s+2*b)/2./(s+b)/TMath::Sqrt(s+b);
 117         Float_t dsigdb=-s/2./(s+b)/TMath::Sqrt(s+b);
 118         errsig=TMath::Sqrt(dsigds*dsigds*errs*errs+dsigdb*dsigdb*errb*errb);
 119       }
 120       fSignificance->SetElement(i,signif);
 121       fErrSignificance->SetElement(i,errsig);
 122     }
 123   }
 124   fSignificance->SetNameTitle("Significance","Significance");
 125   fErrSignificance->SetNameTitle("ErrorOnSignificance","ErrorOnSignificance");
 126 }
 127 //___________________________________________________________________________
 128 AliMultiDimVector* AliSignificanceCalculator::CalculatePurity() const {
 129   // calculates purity
 130   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 131
 132   AliMultiDimVector* purity=new AliMultiDimVector();
 133   purity->CopyStructure(fSignal);
 134   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 135     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 136       Float_t s=fSignal->GetElement(i)*fNormSig;
 137       Float_t b=fBackground->GetElement(i)*fNormBkg;
 138       Float_t pur=0.;
 139       if((s+b)>0) pur=s/(s+b);
 140       purity->SetElement(i,pur);
 141     }
 142   }
 143   purity->SetNameTitle("Purity","Purity");
 144   return purity;
 145 }
 146 //___________________________________________________________________________
 147 AliMultiDimVector* AliSignificanceCalculator::CalculatePurityError() const {
 148   // calculates error on purity
 149   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 150
 151   AliMultiDimVector* epurity=new AliMultiDimVector();
 152   epurity->CopyStructure(fSignal);
 153   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 154     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 155       Float_t s=fSignal->GetElement(i)*fNormSig;
 156       Float_t b=fBackground->GetElement(i)*fNormBkg;
 157       Float_t epur=0.;
 158       if((s+b)>0){
 159         Float_t errs,errb;
 160         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 161         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 162         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 163         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 164         Float_t dpurds=b/(s+b)/(s+b);
 165         Float_t dpurdb=-s/(s+b)/(s+b);
 166         epur=TMath::Sqrt(dpurds*dpurds*errs*errs+dpurdb*dpurdb*errb*errb);
 167       }
 168       epurity->SetElement(i,epur);
 169     }
 170   }
 171   epurity->SetNameTitle("ErrorOnPurity","ErrorOnPurity");
 172   return epurity;
 173 }
 174 //___________________________________________________________________________
 175 AliMultiDimVector* AliSignificanceCalculator::CalculateSOverB() const {
 176   // Signal over Background
 177   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 178
 179   AliMultiDimVector* sob=new AliMultiDimVector();
 180   sob->CopyStructure(fSignal);
 181   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 182     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 183       Float_t s=fSignal->GetElement(i)*fNormSig;
 184       Float_t b=fBackground->GetElement(i)*fNormBkg;
 185       Float_t soverb=0.;
 186       if(b>0) soverb=s/b;
 187       sob->SetElement(i,soverb);
 188     }
 189   }
 190   sob->SetNameTitle("SoverB","SoverB");
 191   return sob;
 192 }
 193 //___________________________________________________________________________
 194 AliMultiDimVector* AliSignificanceCalculator::CalculateSOverBError() const {
 195   // Error on Signal over Background
 196   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 197
 198   AliMultiDimVector* esob=new AliMultiDimVector();
 199   esob->CopyStructure(fSignal);
 200   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 201     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 202       Float_t s=fSignal->GetElement(i)*fNormSig;
 203       Float_t b=fBackground->GetElement(i)*fNormBkg;
 204       Float_t esoverb=0.;
 205       if(b>0){
 206         Float_t soverb=s/b;
 207         Float_t errs,errb;
 208         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 209         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 210         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 211         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 212         esoverb=soverb*TMath::Sqrt(errs*errs/s/s+errb*errb/b/b);
 213       }
 214       esob->SetElement(i,esoverb);
 215     }
 216   }
 217   esob->SetNameTitle("ErrorOnSoverB","ErrorOnSoverB");
 218   return esob;
 219 }