[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliSignificanceCalculator.cxx

/**************************************************************************
 * 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 "AliMultiDimVector.h"
#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;
}
Commit	Line	Data
ac7636a0	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
7b45817b	16	/* $Id$ */
ac7636a0	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
cc4bd7b4	27	#include "AliMultiDimVector.h"
ac7636a0	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(dsigdsdsigdserrserrs+dsigdbdsigdberrberrb);
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(dpurdsdpurdserrserrs+dpurdbdpurdberrberrb);
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=soverbTMath::Sqrt(errserrs/s/s+errb*errb/b/b);
214	}
215	esob->SetElement(i,esoverb);
216	}
217	}
218	esob->SetNameTitle("ErrorOnSoverB","ErrorOnSoverB");
219	return esob;
220	}