[u/mrichter/AliRoot.git] / PWG3 / 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 "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
	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(dsigdsdsigdserrserrs+dsigdbdsigdberrberrb);
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(dpurdsdpurdserrserrs+dpurdbdpurdberrberrb);
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=soverbTMath::Sqrt(errserrs/s/s+errb*errb/b/b);
213	}
214	esob->SetElement(i,esoverb);
215	}
216	}
217	esob->SetNameTitle("ErrorOnSoverB","ErrorOnSoverB");
218	return esob;
219	}