ea361c8963b7e9fc438997b4efa5bac2a04b7ff3
[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[]={0x0};
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::kSPECIESC;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::kSPECIESC;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::kSPECIESC;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::EParticleType)AliPID::kSPECIESC) ) ){
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 (fPriorsDistributions[i]) {
95       delete fPriorsDistributions[i]; 
96     }
97     fPriorsDistributions[i]=new TH1F(*prior);
98   }
99 }
100
101 //-------------------------------------------------------------------------------------------------     
102 UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
103         Double_t pITS[fSelectedSpecies];
104         Double_t pTPC[fSelectedSpecies];
105         Double_t pTRD[fSelectedSpecies];
106         Double_t pTOF[fSelectedSpecies];
107         Double_t pHMPID[fSelectedSpecies];
108         Double_t pEMCAL[fSelectedSpecies];
109         Double_t pPHOS[fSelectedSpecies];
110         for (Int_t i=0;i<fSelectedSpecies;i++) {
111          pITS[i]=1.;
112          pTPC[i]=1.;
113          pTRD[i]=1.;
114          pTOF[i]=1.;
115          pHMPID[i]=1.;
116          pEMCAL[i]=1.;
117          pPHOS[i]=1.;
118         }
119         UInt_t usedMask=0;
120         AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
121         Double_t p[fSelectedSpecies];
122         memset(p,0,sizeof(Double_t)*fSelectedSpecies);
123
124         // getting probability distributions for selected detectors only
125         if (fDetectorMask & AliPIDResponse::kDetITS) {
126           status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
127           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
128         }
129
130         if (fDetectorMask & AliPIDResponse::kDetTPC) { 
131           status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
132           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
133         }
134
135
136         if (fDetectorMask & AliPIDResponse::kDetTRD) { 
137           status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
138           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
139         }
140
141         if (fDetectorMask &  AliPIDResponse::kDetTOF) { 
142           status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
143           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
144         }
145
146         if (fDetectorMask & AliPIDResponse::kDetHMPID) { 
147           status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
148           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
149         }
150
151
152         if (fDetectorMask & AliPIDResponse::kDetEMCAL) { 
153           status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
154           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
155         }
156
157
158         if (fDetectorMask & AliPIDResponse::kDetPHOS) { 
159           status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
160           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
161         }
162
163
164         for (Int_t i =0; i<fSelectedSpecies; i++){
165           p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
166         }
167         Double_t priors[fSelectedSpecies];
168         memset(priors,0,fSelectedSpecies*sizeof(Double_t));
169         if (fEnablePriors){
170           GetPriors(track,priors,response->GetCurrentCentrality());
171           
172           // for the moment we have four cases
173           // TPC+TRD      --> apply TRD propagation factors
174           // TPC+TOF      --> apply TOF propagation factors (TRD may be present)
175           // TPC+EMCAL    --> apply EMCAL propagation factors (TRD and TOF if present)
176           if(fUseDefaultTPCPriors) {
177              Double_t pt=TMath::Abs(track->Pt());
178              if ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL ) { // EMCAL is the outer having prop. factors for the moment
179                // EMCal case (for the moment only in combination with TPC)
180                // propagation factors determined from LHC11d MC (LHC12a15f)
181                // v2 clusterizer, dEta < 0.015, dPhi < 0.03, NonLinearityFunction = 6
182                
183                Float_t electronEMCALfactor=0.1;
184                Float_t kaonEMCALfactor = 0.1;
185                Float_t protonEMCALfactor = 0.1;
186                
187                if(track->Charge() > 0){
188                  // positiv charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
189                  if(pt > 0.75 && pt < 20.0){
190                    electronEMCALfactor =  ( 0.214159 * ( 1 - TMath::Exp(-TMath::Power(pt,0.484512)/0.700499)*TMath::Power(pt,-0.669644)) ) /  ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
191                    kaonEMCALfactor =  ( 0.208686 * ( 1 - TMath::Exp(-TMath::Power(pt,-3.98149e-05)/1.28447)*TMath::Power(pt,-0.629191)) ) /  ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
192                    protonEMCALfactor =  ( 0.27555 * ( 1 - TMath::Exp(-TMath::Power(pt,-1.97226e-05)/1.52719)*TMath::Power(pt,-0.209068)) ) /  ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
193                    
194                    }
195                }
196                
197                if(track->Charge() < 0){
198                  // negative charge  (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
199                  if(pt > 0.75 && pt < 20.0){               
200                    electronEMCALfactor =  ( 0.216895 * ( 1 - TMath::Exp(-TMath::Power(pt,0.000105924)/0.865938)*TMath::Power(pt,-1.32787)) ) /  ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
201                    kaonEMCALfactor =  ( 0.204117 * ( 1 - TMath::Exp(-TMath::Power(pt,-1.6853e-05)/1.61765)*TMath::Power(pt,-0.738355)) ) /  ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
202                    protonEMCALfactor =  ( 0.215679 * ( 1 - TMath::Exp(-TMath::Power(pt,-4.10015e-05)/1.40921)*TMath::Power(pt,-0.533752)) ) /  ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
203                  }
204                  }
205                p[0] *= electronEMCALfactor;
206                p[3] *= kaonEMCALfactor;
207                p[4] *= protonEMCALfactor;             
208              } // end of EMCAL case
209              else if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){
210                Float_t kaonTOFfactor = 0.1;
211                if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
212                Float_t protonTOFfactor = 0.1;
213                if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
214                
215                if(track->Charge() < 0){ // for negative tracks
216                  kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
217                  protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
218                }
219                // TODO: we may need an electron factor for TOF as well, especially if TRD is present!
220                p[3] *= kaonTOFfactor;
221                p[4] *= protonTOFfactor;
222              }  // end of TOF case
223              else if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
224                  Float_t electronTRDfactor=0.1;
225                  Float_t kaonTRDfactor = 0.1;
226                  Float_t protonTRDfactor = 0.1;
227                  
228                  if(track->Charge() > 0){
229                    // positiv charge
230                    if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
231                    if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
232                    if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
233                  }
234                  
235                  if(track->Charge() < 0){
236                    // negative charge
237                    if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
238                    if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
239                    if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
240                  }
241                  // what about electrons
242                  p[0] *= electronTRDfactor;
243                  p[3] *= kaonTRDfactor;
244                  p[4] *= protonTRDfactor;             
245              } // end of TRD case
246           } // end of fUseDefaultTPCPriors
247         }   // end of use priors
248         else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
249         ComputeBayesProbabilities(bayesProbabilities,p,priors);
250         return usedMask;
251 }
252
253
254 //-------------------------------------------------------------------------------------------------
255 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
256         
257         //
258         // get priors from distributions
259         //
260         
261         Double_t pt=TMath::Abs(track->Pt());
262
263         if(fUseDefaultTPCPriors){ // use default priors if requested
264           Float_t usedCentr = centrality+5*(centrality>0);
265           if(usedCentr < -0.99) usedCentr = -0.99;
266           else if(usedCentr > 99) usedCentr = 99;
267           if(pt > 9.99) pt = 9.99;
268           else if(pt < 0.01)  pt = 0.01;
269           
270           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
271
272           return;
273         }
274         
275         Double_t sumPriors = 0;
276         for (Int_t i=0;i<fSelectedSpecies;++i){
277                 if (fPriorsDistributions[i]){
278                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
279                 }
280                 else {
281                         p[i]=0.;
282                 }
283                 sumPriors+=p[i];                
284         }
285
286         // normalizing priors
287         if (sumPriors == 0) return;
288         for (Int_t i=0;i<fSelectedSpecies;++i){
289            p[i]=p[i]/sumPriors;
290         }
291         return;
292 }
293 //-------------------------------------------------------------------------------------------------
294 void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
295         
296         //
297         // get priors from distributions
298         //
299         Double_t pt=TMath::Abs(track->Pt());
300         
301         if(fUseDefaultTPCPriors){ // use default priors if requested
302           Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
303           if(usedCentr < -0.99) usedCentr = -0.99;
304           else if(usedCentr > 99) usedCentr = 99;
305           if(pt > 9.99) pt = 9.99;
306           else if(pt < 0.01)  pt = 0.01;
307           
308           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
309
310           // Extra factor if TOF matching was required
311           if(detUsed & AliPIDResponse::kDetTOF){
312             Float_t kaonTOFfactor = 0.1;
313             if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
314             Float_t protonTOFfactor = 0.1;
315             if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
316             
317             if(track->Charge() < 0){ // for negative tracks
318               kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
319               protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
320             }
321             
322             p[3] *= kaonTOFfactor;
323             p[4] *= protonTOFfactor;
324           }
325
326           return;
327         }
328
329
330         Double_t sumPriors = 0;
331         for (Int_t i=0;i<fSelectedSpecies;++i){
332                 if (fPriorsDistributions[i]){
333                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
334                 }
335                 else {
336                         p[i]=0.;
337                 }
338                 sumPriors+=p[i];                
339         }
340
341         // normalizing priors
342         if (sumPriors == 0) return;
343         for (Int_t i=0;i<fSelectedSpecies;++i){
344            p[i]=p[i]/sumPriors;
345         }
346         return;
347 }
348 //-------------------------------------------------------------------------------------------------     
349 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
350
351
352   //
353   // calculate Bayesian probabilities
354   //
355   Double_t sum = 0.;
356   for (Int_t i = 0; i < fSelectedSpecies; i++) {
357     sum += probDensity[i] * prior[i];
358   }
359   if (sum <= 0) {
360
361     AliError("Invalid probability densities or priors");
362     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
363     return;
364   }
365   for (Int_t i = 0; i < fSelectedSpecies; i++) {
366     probabilities[i] = probDensity[i] * prior[i] / sum;
367   }
368
369
370 }
371
372 //----------------------------------------------------------------------------------------
373 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
374   switch (status) {
375     case AliPIDResponse::kDetNoParams:
376     case AliPIDResponse::kDetNoSignal:
377     break;
378   case AliPIDResponse::kDetPidOk:
379     *maskDetIn = *maskDetIn | bit;
380     break;
381   case AliPIDResponse::kDetMismatch:
382     if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
383     else *maskDetIn = *maskDetIn | bit;
384     break;
385   }
386
387 }
388
389 //----------------------------------------------------------------------------------------
390 void AliPIDCombined::SetDefaultTPCPriors(){
391   fEnablePriors=kTRUE;
392   fUseDefaultTPCPriors = kTRUE;
393
394   //check if priors are already initialized
395   if (fDefaultPriorsTPC[0]) return;
396   
397   TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
398   oadbfilename += "/PIDdefaultPriors.root";
399   TFile * foadb = TFile::Open(oadbfilename.Data());
400   if(!foadb || !foadb->IsOpen()) {
401     delete foadb;
402     AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
403     return;
404   }
405   
406   AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
407   if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
408   
409   const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
410   char name[100];
411
412   for(Int_t i=0;i < AliPID::kSPECIES;i++){
413     snprintf(name,99,"hDefault%sPriors",namespecies[i]);
414     fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
415     if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
416   }
417
418   delete foadb;
419 }