Fixes for #89479: change in AliPIDResponse and AliPIDCombined (STEERBase). F.Noferini
[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 #include "TMath.h"
36 #include "TFile.h"
37
38 #include "AliOADBContainer.h"
39
40 // initialize static members
41 TH2F* AliPIDCombined::fDefaultPriorsTPC[5];
42
43 ClassImp(AliPIDCombined);
44
45 AliPIDCombined::AliPIDCombined():
46         TNamed("CombinedPID","CombinedPID"),
47         fDetectorMask(0),
48         fRejectMismatchMask(0x7F),
49         fEnablePriors(kTRUE),
50         fSelectedSpecies(AliPID::kSPECIES),
51         fUseDefaultTPCPriors(kFALSE)
52 {       
53   //
54   // default constructor
55   //    
56   for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
57   AliLog::SetClassDebugLevel("AliPIDCombined",10);
58 }
59
60 //-------------------------------------------------------------------------------------------------     
61 AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
62         TNamed(name,title),
63         fDetectorMask(0),
64         fRejectMismatchMask(0x7F),
65         fEnablePriors(kTRUE),
66         fSelectedSpecies(AliPID::kSPECIES),
67         fUseDefaultTPCPriors(kFALSE)
68 {
69   //
70   // default constructor with name and title
71   //
72   for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
73   AliLog::SetClassDebugLevel("AliPIDCombined",10);
74
75 }
76
77 //-------------------------------------------------------------------------------------------------     
78 AliPIDCombined::~AliPIDCombined() {
79
80   for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
81     if(fPriorsDistributions[i])
82       delete fPriorsDistributions[i];
83   }
84 }
85
86 //-------------------------------------------------------------------------------------------------     
87 void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
88   if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
89     AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
90     return;
91   }
92   if(prior) {
93     Int_t i = (Int_t)type;
94     if (i >= AliPID::kSPECIES ) {
95       if (i < (Int_t)AliPID::kDeuteron) {
96         AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
97         return;
98       } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
99     }
100     if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) {             // coverity fix (to put it mildly....)
101       AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
102       return;
103     }
104     if (fPriorsDistributions[i]) {
105       delete fPriorsDistributions[i]; 
106     }
107     fPriorsDistributions[i]=new TH1F(*prior);
108   }
109 }
110
111 //-------------------------------------------------------------------------------------------------     
112 UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
113         Double_t pITS[fSelectedSpecies];
114         Double_t pTPC[fSelectedSpecies];
115         Double_t pTRD[fSelectedSpecies];
116         Double_t pTOF[fSelectedSpecies];
117         Double_t pHMPID[fSelectedSpecies];
118         Double_t pEMCAL[fSelectedSpecies];
119         Double_t pPHOS[fSelectedSpecies];
120         for (Int_t i=0;i<fSelectedSpecies;i++) {
121          pITS[i]=1.;
122          pTPC[i]=1.;
123          pTRD[i]=1.;
124          pTOF[i]=1.;
125          pHMPID[i]=1.;
126          pEMCAL[i]=1.;
127          pPHOS[i]=1.;
128         }
129         UInt_t usedMask=0;
130         AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
131         Double_t p[fSelectedSpecies];
132
133         // getting probability distributions for selected detectors only
134         if (fDetectorMask & AliPIDResponse::kDetITS) {
135           status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
136           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
137         }
138
139         if (fDetectorMask & AliPIDResponse::kDetTPC) { 
140           status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
141           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
142         }
143
144
145         if (fDetectorMask & AliPIDResponse::kDetTRD) { 
146           status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
147           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
148         }
149
150         if (fDetectorMask &  AliPIDResponse::kDetTOF) { 
151           status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
152           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
153         }
154
155         if (fDetectorMask & AliPIDResponse::kDetHMPID) { 
156           status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
157           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
158         }
159
160
161         if (fDetectorMask & AliPIDResponse::kDetEMCAL) { 
162           status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
163           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
164         }
165
166
167         if (fDetectorMask & AliPIDResponse::kDetPHOS) { 
168           status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
169           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
170         }
171
172
173         for (Int_t i =0; i<fSelectedSpecies; i++){
174           p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
175         }
176         Double_t priors[fSelectedSpecies];
177         if (fEnablePriors){
178           GetPriors(track,priors,response->GetCurrentCentrality());
179           
180           // apply tof matching efficiency to priors if TOF joined PID for this track
181           if(fUseDefaultTPCPriors && (usedMask & AliPIDResponse::kDetTOF)){
182             Double_t pt=TMath::Abs(track->Pt());
183             Float_t kaonTOFfactor = 0.1;
184             if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
185             Float_t protonTOFfactor = 0.1;
186             if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
187             
188             if(track->Charge() < 0){ // for negative tracks
189               kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
190               protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
191             }
192             
193             p[3] *= kaonTOFfactor;
194             p[4] *= protonTOFfactor;
195           }
196         }
197         else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
198         ComputeBayesProbabilities(bayesProbabilities,p,priors);
199         return usedMask;
200 }
201
202
203 //-------------------------------------------------------------------------------------------------
204 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
205         
206         //
207         // get priors from distributions
208         //
209         
210         Double_t pt=TMath::Abs(track->Pt());
211
212         if(fUseDefaultTPCPriors){ // use default priors if requested
213           Float_t usedCentr = centrality+5*(centrality>0);
214           if(usedCentr < -0.99) usedCentr = -0.99;
215           else if(usedCentr > 99) usedCentr = 99;
216           if(pt > 9.99) pt = 9.99;
217           else if(pt < 0.01)  pt = 0.01;
218           
219           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
220
221           return;
222         }
223         
224         Double_t sumPriors = 0;
225         for (Int_t i=0;i<fSelectedSpecies;++i){
226                 if (fPriorsDistributions[i]){
227                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
228                 }
229                 else {
230                         p[i]=0.;
231                 }
232                 sumPriors+=p[i];                
233         }
234
235         // normalizing priors
236         if (sumPriors == 0) return;
237         for (Int_t i=0;i<fSelectedSpecies;++i){
238            p[i]=p[i]/sumPriors;
239         }
240         return;
241 }
242 //-------------------------------------------------------------------------------------------------
243 void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
244         
245         //
246         // get priors from distributions
247         //
248         Double_t pt=TMath::Abs(track->Pt());
249         
250         if(fUseDefaultTPCPriors){ // use default priors if requested
251           Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
252           if(usedCentr < -0.99) usedCentr = -0.99;
253           else if(usedCentr > 99) usedCentr = 99;
254           if(pt > 9.99) pt = 9.99;
255           else if(pt < 0.01)  pt = 0.01;
256           
257           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
258
259           // Extra factor if TOF matching was required
260           if(detUsed & AliPIDResponse::kDetTOF){
261             Float_t kaonTOFfactor = 0.1;
262             if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
263             Float_t protonTOFfactor = 0.1;
264             if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
265             
266             if(track->Charge() < 0){ // for negative tracks
267               kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
268               protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
269             }
270             
271             p[3] *= kaonTOFfactor;
272             p[4] *= protonTOFfactor;
273           }
274
275           return;
276         }
277
278
279         Double_t sumPriors = 0;
280         for (Int_t i=0;i<fSelectedSpecies;++i){
281                 if (fPriorsDistributions[i]){
282                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
283                 }
284                 else {
285                         p[i]=0.;
286                 }
287                 sumPriors+=p[i];                
288         }
289
290         // normalizing priors
291         if (sumPriors == 0) return;
292         for (Int_t i=0;i<fSelectedSpecies;++i){
293            p[i]=p[i]/sumPriors;
294         }
295         return;
296 }
297 //-------------------------------------------------------------------------------------------------     
298 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
299
300
301   //
302   // calculate Bayesian probabilities
303   //
304   Double_t sum = 0.;
305   for (Int_t i = 0; i < fSelectedSpecies; i++) {
306     sum += probDensity[i] * prior[i];
307   }
308   if (sum <= 0) {
309
310     AliError("Invalid probability densities or priors");
311     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
312     return;
313   }
314   for (Int_t i = 0; i < fSelectedSpecies; i++) {
315     probabilities[i] = probDensity[i] * prior[i] / sum;
316   }
317
318
319 }
320
321 //----------------------------------------------------------------------------------------
322 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
323   switch (status) {
324   case AliPIDResponse::kDetNoSignal:
325     break;
326   case AliPIDResponse::kDetPidOk:
327     *maskDetIn = *maskDetIn | bit;
328     break;
329   case AliPIDResponse::kDetMismatch:
330     if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
331     else *maskDetIn = *maskDetIn | bit;
332     break;
333   }
334
335 }
336
337 //----------------------------------------------------------------------------------------
338 void AliPIDCombined::SetDefaultTPCPriors(){
339   fEnablePriors=kTRUE;
340   fUseDefaultTPCPriors = kTRUE;
341
342   TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
343   oadbfilename += "/PIDdefaultPriors.root";
344   TFile * foadb = TFile::Open(oadbfilename.Data());
345   if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
346   
347   AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
348   if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
349   
350   const char *namespecies[5] = {"El","Mu","Pi","Ka","Pr"};
351   char name[100];
352
353   for(Int_t i=0;i < 5;i++){
354     sprintf(name,"hDefault%sPriors",namespecies[i]);
355     fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
356     if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
357   }
358 }