[u/mrichter/AliRoot.git] / STEER / AliPIDCombined.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/


//-----------------------------------------------------------------
//        Base class for combining PID of different detectors    //
//        (user selected) and compute Bayesian probabilities     //
//                                                               //
//                                                               //
//   Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch  //
//                                                               //
//-----------------------------------------------------------------

#include <TH1.h>

#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliPIDResponse.h>

#include "AliPIDCombined.h"

AliPIDCombined::AliPIDCombined():
	TNamed("CombinedPID","CombinedPID"),
	fDetectorMask(0),
	fRejectMismatchMask(0x7F),
	fEnablePriors(kTRUE),
	fSelectedSpecies(AliPID::kSPECIES)
{	
  //
  // default constructor
  //	
  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
  AliLog::SetClassDebugLevel("AliPIDCombined",10);
}

//-------------------------------------------------------------------------------------------------	
AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
	TNamed(name,title),
	fDetectorMask(0),
	fRejectMismatchMask(0x7F),
	fEnablePriors(kTRUE),
	fSelectedSpecies(AliPID::kSPECIES)
{
  //
  // default constructor with name and title
  //
  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
  AliLog::SetClassDebugLevel("AliPIDCombined",10);

}

//-------------------------------------------------------------------------------------------------	
AliPIDCombined::~AliPIDCombined() {

  for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
    if(fPriorsDistributions[i])
      delete fPriorsDistributions[i];
  }
}

//-------------------------------------------------------------------------------------------------	
void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
  if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
    AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
    return;
  }
  if(prior) {
    Int_t i = (Int_t)type;
    if (i >= AliPID::kSPECIES ) {
      if (i < (Int_t)AliPID::kDeuteron) {
	AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
	return;
      } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
    }
    if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) {             // coverity fix (to put it mildly....)
      AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
      return;
    }
    if (fPriorsDistributions[i]) {
      delete fPriorsDistributions[i]; 
    }
    fPriorsDistributions[i]=new TH1F(*prior);
  }
}

//-------------------------------------------------------------------------------------------------	
UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
	Double_t pITS[fSelectedSpecies];
	Double_t pTPC[fSelectedSpecies];
	Double_t pTRD[fSelectedSpecies];
	Double_t pTOF[fSelectedSpecies];
	Double_t pHMPID[fSelectedSpecies];
	Double_t pEMCAL[fSelectedSpecies];
	Double_t pPHOS[fSelectedSpecies];
	for (Int_t i=0;i<fSelectedSpecies;i++) {
	 pITS[i]=1.;
	 pTPC[i]=1.;
	 pTRD[i]=1.;
	 pTOF[i]=1.;
	 pHMPID[i]=1.;
	 pEMCAL[i]=1.;
	 pPHOS[i]=1.;
	}
	UInt_t usedMask=0;
	AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
	Double_t p[fSelectedSpecies];

	// getting probability distributions for selected detectors only
	if (fDetectorMask & AliPIDResponse::kDetITS) {
	  status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
	}

	if (fDetectorMask & AliPIDResponse::kDetTPC) { 
	  status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
	}


	if (fDetectorMask & AliPIDResponse::kDetTRD) { 
	  status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
	}

	if (fDetectorMask &  AliPIDResponse::kDetTOF) { 
	  status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
	}

	if (fDetectorMask & AliPIDResponse::kDetHMPID) { 
	  status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
	}


	if (fDetectorMask & AliPIDResponse::kDetEMCAL) { 
	  status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
	}


	if (fDetectorMask & AliPIDResponse::kDetPHOS) { 
	  status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
	}


	for (Int_t i =0; i<fSelectedSpecies; i++){
	  p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
	}
	Double_t priors[fSelectedSpecies];
	if (fEnablePriors) GetPriors(track,priors);
	else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
	ComputeBayesProbabilities(bayesProbabilities,p,priors);
	return usedMask;
}


//-------------------------------------------------------------------------------------------------
void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {
	
	//
	// get priors from distributions
	//
	
	Double_t pt=TMath::Abs(track->Pt());
	Double_t sumPriors = 0;
	for (Int_t i=0;i<fSelectedSpecies;++i){
		if (fPriorsDistributions[i]){
			p[i]=fPriorsDistributions[i]->Interpolate(pt);
		}
		else {
			p[i]=0.;
		}
		sumPriors+=p[i];		
	}

	// normalizing priors
	if (sumPriors == 0) return;
	for (Int_t i=0;i<fSelectedSpecies;++i){
	   p[i]=p[i]/sumPriors;
	}
	return;
}

//-------------------------------------------------------------------------------------------------	
void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {


  //
  // calculate Bayesian probabilities
  //
  Double_t sum = 0.;
  for (Int_t i = 0; i < fSelectedSpecies; i++) {
    sum += probDensity[i] * prior[i];
  }
  if (sum <= 0) {
    AliError("Invalid probability densities or priors");
    for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
    return;
  }
  for (Int_t i = 0; i < fSelectedSpecies; i++) {
    probabilities[i] = probDensity[i] * prior[i] / sum;
  }


}

//----------------------------------------------------------------------------------------
void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
  switch (status) {
  case AliPIDResponse::kDetNoSignal:
    break;
  case AliPIDResponse::kDetPidOk:
    *maskDetIn = *maskDetIn | bit;
    break;
  case AliPIDResponse::kDetMismatch:
    if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
    else *maskDetIn = *maskDetIn | bit;
    break;
  }

}

ClassImp(AliPIDCombined);
Commit	Line	Data
b8bfee30	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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
	17	//-----------------------------------------------------------------
	18	// Base class for combining PID of different detectors //
	19	// (user selected) and compute Bayesian probabilities //
	20	// //
	21	// //
	22	// Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch //
	23	// //
	24	//-----------------------------------------------------------------
	25
	26	#include <TH1.h>
	27
	28	#include <AliVTrack.h>
	29	#include <AliLog.h>
	30	#include <AliPID.h>
	31	#include <AliPIDResponse.h>
	32
	33	#include "AliPIDCombined.h"
	34
	35	AliPIDCombined::AliPIDCombined():
	36	TNamed("CombinedPID","CombinedPID"),
	37	fDetectorMask(0),
	38	fRejectMismatchMask(0x7F),
	39	fEnablePriors(kTRUE),
	40	fSelectedSpecies(AliPID::kSPECIES)
	41	{
	42	//
	43	// default constructor
	44	//
	45	for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
	46	AliLog::SetClassDebugLevel("AliPIDCombined",10);
	47	}
	48
	49	//-------------------------------------------------------------------------------------------------
	50	AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
	51	TNamed(name,title),
	52	fDetectorMask(0),
	53	fRejectMismatchMask(0x7F),
	54	fEnablePriors(kTRUE),
	55	fSelectedSpecies(AliPID::kSPECIES)
	56	{
	57	//
	58	// default constructor with name and title
	59	//
	60	for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
	61	AliLog::SetClassDebugLevel("AliPIDCombined",10);
	62
	63	}
	64
65	//-------------------------------------------------------------------------------------------------
66	AliPIDCombined::~AliPIDCombined() {
67
68	for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
69	if(fPriorsDistributions[i])
70	delete fPriorsDistributions[i];
71	}
72	}
73
74	//-------------------------------------------------------------------------------------------------
75	void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
76	if ( (type < 0) \|\| ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
77	AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
78	return;
79	}
80	if(prior) {
8de3b790	81	Int_t i = (Int_t)type;
	82	if (i >= AliPID::kSPECIES ) {
	83	if (i < (Int_t)AliPID::kDeuteron) {
	84	AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
b8bfee30	85	return;
8de3b790	86	} else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
	87	}
	88	if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) { // coverity fix (to put it mildly....)
	89	AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
	90	return;
b8bfee30	91	}
	92	if (fPriorsDistributions[i]) {
	93	delete fPriorsDistributions[i];
	94	}
	95	fPriorsDistributions[i]=new TH1F(*prior);
	96	}
	97	}
	98
	99	//-------------------------------------------------------------------------------------------------
	100	UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack track, const AliPIDResponse response, Double_t* bayesProbabilities) const {
	101	Double_t pITS[fSelectedSpecies];
	102	Double_t pTPC[fSelectedSpecies];
	103	Double_t pTRD[fSelectedSpecies];
	104	Double_t pTOF[fSelectedSpecies];
	105	Double_t pHMPID[fSelectedSpecies];
	106	Double_t pEMCAL[fSelectedSpecies];
	107	Double_t pPHOS[fSelectedSpecies];
	108	for (Int_t i=0;i<fSelectedSpecies;i++) {
	109	pITS[i]=1.;
	110	pTPC[i]=1.;
	111	pTRD[i]=1.;
	112	pTOF[i]=1.;
	113	pHMPID[i]=1.;
	114	pEMCAL[i]=1.;
	115	pPHOS[i]=1.;
	116	}
	117	UInt_t usedMask=0;
	118	AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
	119	Double_t p[fSelectedSpecies];
	120
	121	// getting probability distributions for selected detectors only
	122	if (fDetectorMask & AliPIDResponse::kDetITS) {
	123	status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
	124	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
	125	}
	126
	127	if (fDetectorMask & AliPIDResponse::kDetTPC) {
	128	status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
	129	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
	130	}
	131
	132
	133	if (fDetectorMask & AliPIDResponse::kDetTRD) {
	134	status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
	135	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
	136	}
	137
	138	if (fDetectorMask & AliPIDResponse::kDetTOF) {
	139	status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
	140	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
	141	}
	142
	143	if (fDetectorMask & AliPIDResponse::kDetHMPID) {
	144	status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
	145	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
	146	}
	147
	148
	149	if (fDetectorMask & AliPIDResponse::kDetEMCAL) {
	150	status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
	151	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
	152	}
	153
	154
155	if (fDetectorMask & AliPIDResponse::kDetPHOS) {
156	status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
157	SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
158	}
159
160
161	for (Int_t i =0; i<fSelectedSpecies; i++){
162	p[i] = pITS[i]pTPC[i]pTRD[i]pTOF[i]pHMPID[i]pEMCAL[i]pPHOS[i];
163	}
164	Double_t priors[fSelectedSpecies];
165	if (fEnablePriors) GetPriors(track,priors);
166	else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
167	ComputeBayesProbabilities(bayesProbabilities,p,priors);
168	return usedMask;
169	}
170
171
172	//-------------------------------------------------------------------------------------------------
173	void AliPIDCombined::GetPriors(const AliVTrack track, Double_t p) const {
174
175	//
176	// get priors from distributions
177	//
178
179	Double_t pt=TMath::Abs(track->Pt());
180	Double_t sumPriors = 0;
181	for (Int_t i=0;i<fSelectedSpecies;++i){
182	if (fPriorsDistributions[i]){
183	p[i]=fPriorsDistributions[i]->Interpolate(pt);
184	}
185	else {
186	p[i]=0.;
187	}
188	sumPriors+=p[i];
189	}
190
191	// normalizing priors
192	if (sumPriors == 0) return;
193	for (Int_t i=0;i<fSelectedSpecies;++i){
194	p[i]=p[i]/sumPriors;
195	}
196	return;
197	}
198
199	//-------------------------------------------------------------------------------------------------
200	void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
201
202
203	//
204	// calculate Bayesian probabilities
205	//
206	Double_t sum = 0.;
207	for (Int_t i = 0; i < fSelectedSpecies; i++) {
208	sum += probDensity[i] * prior[i];
209	}
210	if (sum <= 0) {
211	AliError("Invalid probability densities or priors");
212	for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
213	return;
214	}
215	for (Int_t i = 0; i < fSelectedSpecies; i++) {
216	probabilities[i] = probDensity[i] * prior[i] / sum;
217	}
218
219
220	}
221
222	//----------------------------------------------------------------------------------------
223	void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t maskDetIn, AliPIDResponse::EDetCode bit, Double_t p) const {
224	switch (status) {
225	case AliPIDResponse::kDetNoSignal:
226	break;
227	case AliPIDResponse::kDetPidOk:
228	maskDetIn = maskDetIn \| bit;
229	break;
230	case AliPIDResponse::kDetMismatch:
231	if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
232	else maskDetIn = maskDetIn \| bit;
233	break;
234	}
235
236	}
237
238	ClassImp(AliPIDCombined);
239