PWGHF/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 "AliMultiDimVector.h"
  28 #include "AliSignificanceCalculator.h"
  29 #include "TMath.h"
  30 #include "AliLog.h"
  31
  32 ClassImp(AliSignificanceCalculator)
  33 //___________________________________________________________________________
  34 AliSignificanceCalculator::AliSignificanceCalculator():
  35 fSignal(0),
  36 fErrSquareSignal(0),
  37 fBackground(0),
  38 fErrSquareBackground(0),
  39 fSignificance(0),
  40 fErrSignificance(0),
  41 fNormSig(1.),
  42 fNormBkg(1.)
  43 {
  44   // default constructor
  45   fSignal=new AliMultiDimVector();
  46   fBackground=new AliMultiDimVector();
  47 }
  48 //___________________________________________________________________________
  49 AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, Float_t normsig, Float_t normbkg):
  50 fSignal(sig),
  51 fErrSquareSignal(0),
  52 fBackground(bkg),
  53 fErrSquareBackground(0),
  54 fSignificance(0),
  55 fErrSignificance(0),
  56 fNormSig(normsig),
  57 fNormBkg(normbkg)
  58 {
  59   // standard constructor
  60   if(fSignal && fBackground) CalculateSignificance();
  61 }
  62 //___________________________________________________________________________
  63 AliSignificanceCalculator::AliSignificanceCalculator(AliMultiDimVector* sig, AliMultiDimVector* bkg, AliMultiDimVector* err2sig, AliMultiDimVector* err2bkg, Float_t normsig, Float_t normbkg):
  64 fSignal(sig),
  65 fErrSquareSignal(err2sig),
  66 fBackground(bkg),
  67 fErrSquareBackground(err2bkg),
  68 fSignificance(0),
  69 fErrSignificance(0),
  70 fNormSig(normsig),
  71 fNormBkg(normbkg)
  72 {
  73   // standard constructor
  74   if(fSignal && fBackground) CalculateSignificance();
  75 }
  76 //___________________________________________________________________________
  77 AliSignificanceCalculator::~AliSignificanceCalculator(){
  78   // destructor
  79   if(fSignal) delete fSignal;
  80   if(fBackground) delete fBackground;
  81   if(fSignificance) delete fSignificance;
  82   if(fErrSignificance) delete fErrSignificance;
  83 }
  84 //___________________________________________________________________________
  85 Bool_t AliSignificanceCalculator::Check() const {
  86   // checks AliMultiDimVector dimension and normalization
  87   if(fSignal==0 || fBackground==0) return kFALSE;
  88   if(fNormSig==0. || fNormBkg==0.) return kFALSE;
  89   if(fSignal->GetNTotCells() != fBackground->GetNTotCells()) return kFALSE;
  90   return kTRUE;
  91 }
  92 //___________________________________________________________________________
  93 void AliSignificanceCalculator::CalculateSignificance(){
  94   // calculates significance and its error
  95   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
  96
  97   if(fSignificance) delete fSignificance;
  98   if(fErrSignificance) delete fErrSignificance;
  99   fSignificance=new AliMultiDimVector();
 100   fSignificance->CopyStructure(fSignal);
 101   fErrSignificance=new AliMultiDimVector();
 102   fErrSignificance->CopyStructure(fSignal);
 103
 104   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 105     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 106       Float_t s=fSignal->GetElement(i)*fNormSig;
 107       Float_t b=fBackground->GetElement(i)*fNormBkg;
 108       Float_t signif=0.;
 109       Float_t errsig=0.;
 110       if((s+b)>0){
 111         signif=s/TMath::Sqrt(s+b);
 112         Float_t errs,errb;
 113         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 114         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 115         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 116         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 117         Float_t dsigds=(s+2*b)/2./(s+b)/TMath::Sqrt(s+b);
 118         Float_t dsigdb=-s/2./(s+b)/TMath::Sqrt(s+b);
 119         errsig=TMath::Sqrt(dsigds*dsigds*errs*errs+dsigdb*dsigdb*errb*errb);
 120       }
 121       fSignificance->SetElement(i,signif);
 122       fErrSignificance->SetElement(i,errsig);
 123     }
 124   }
 125   fSignificance->SetNameTitle("Significance","Significance");
 126   fErrSignificance->SetNameTitle("ErrorOnSignificance","ErrorOnSignificance");
 127 }
 128 //___________________________________________________________________________
 129 AliMultiDimVector* AliSignificanceCalculator::CalculatePurity() const {
 130   // calculates purity
 131   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 132
 133   AliMultiDimVector* purity=new AliMultiDimVector();
 134   purity->CopyStructure(fSignal);
 135   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 136     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 137       Float_t s=fSignal->GetElement(i)*fNormSig;
 138       Float_t b=fBackground->GetElement(i)*fNormBkg;
 139       Float_t pur=0.;
 140       if((s+b)>0) pur=s/(s+b);
 141       purity->SetElement(i,pur);
 142     }
 143   }
 144   purity->SetNameTitle("Purity","Purity");
 145   return purity;
 146 }
 147 //___________________________________________________________________________
 148 AliMultiDimVector* AliSignificanceCalculator::CalculatePurityError() const {
 149   // calculates error on purity
 150   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 151
 152   AliMultiDimVector* epurity=new AliMultiDimVector();
 153   epurity->CopyStructure(fSignal);
 154   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 155     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 156       Float_t s=fSignal->GetElement(i)*fNormSig;
 157       Float_t b=fBackground->GetElement(i)*fNormBkg;
 158       Float_t epur=0.;
 159       if((s+b)>0){
 160         Float_t errs,errb;
 161         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 162         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 163         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 164         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 165         Float_t dpurds=b/(s+b)/(s+b);
 166         Float_t dpurdb=-s/(s+b)/(s+b);
 167         epur=TMath::Sqrt(dpurds*dpurds*errs*errs+dpurdb*dpurdb*errb*errb);
 168       }
 169       epurity->SetElement(i,epur);
 170     }
 171   }
 172   epurity->SetNameTitle("ErrorOnPurity","ErrorOnPurity");
 173   return epurity;
 174 }
 175 //___________________________________________________________________________
 176 AliMultiDimVector* AliSignificanceCalculator::CalculateSOverB() const {
 177   // Signal over Background
 178   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 179
 180   AliMultiDimVector* sob=new AliMultiDimVector();
 181   sob->CopyStructure(fSignal);
 182   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 183     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 184       Float_t s=fSignal->GetElement(i)*fNormSig;
 185       Float_t b=fBackground->GetElement(i)*fNormBkg;
 186       Float_t soverb=0.;
 187       if(b>0) soverb=s/b;
 188       sob->SetElement(i,soverb);
 189     }
 190   }
 191   sob->SetNameTitle("SoverB","SoverB");
 192   return sob;
 193 }
 194 //___________________________________________________________________________
 195 AliMultiDimVector* AliSignificanceCalculator::CalculateSOverBError() const {
 196   // Error on Signal over Background
 197   if(!Check()) AliFatal("Signal and Background AliMultiDimVector dimensions do not match!");
 198
 199   AliMultiDimVector* esob=new AliMultiDimVector();
 200   esob->CopyStructure(fSignal);
 201   for(ULong64_t i=0;i<fSignal->GetNTotCells();i++) {
 202     if(fSignal->GetElement(i)!=-1 && fBackground->GetElement(i)!=-1){
 203       Float_t s=fSignal->GetElement(i)*fNormSig;
 204       Float_t b=fBackground->GetElement(i)*fNormBkg;
 205       Float_t esoverb=0.;
 206       if(b>0){
 207         Float_t soverb=s/b;
 208         Float_t errs,errb;
 209         if(fErrSquareSignal) errs=TMath::Sqrt(fErrSquareSignal->GetElement(i))*fNormSig;
 210         else errs=TMath::Sqrt(fSignal->GetElement(i))*fNormSig; // Poisson statistics
 211         if(fErrSquareBackground) errb=TMath::Sqrt(fErrSquareBackground->GetElement(i))*fNormBkg;
 212         else errb=TMath::Sqrt(fBackground->GetElement(i))*fNormBkg; // Poisson
 213         esoverb=soverb*TMath::Sqrt(errs*errs/s/s+errb*errb/b/b);
 214       }
 215       esob->SetElement(i,esoverb);
 216     }
 217   }
 218   esob->SetNameTitle("ErrorOnSoverB","ErrorOnSoverB");
 219   return esob;
 220 }