Updates in PbPb cuts (Andrea Rossi)
[u/mrichter/AliRoot.git] / 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 }