Moving the classes that belong to the following libraries: STEERBase, ESD, CDB, AOD...
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDCombined.cxx
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) {
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));
85         return;
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;
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