Line data Source code
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 : TH2F* AliPIDCombined::fDefaultPriorsTPCpPb[]={0x0};
43 : Float_t AliPIDCombined::fTOFmismProb = 0;
44 :
45 176 : ClassImp(AliPIDCombined);
46 :
47 : AliPIDCombined::AliPIDCombined():
48 0 : TNamed("CombinedPID","CombinedPID"),
49 0 : fDetectorMask(0),
50 0 : fRejectMismatchMask(0x7F),
51 0 : fEnablePriors(kFALSE),
52 0 : fSelectedSpecies(AliPID::kSPECIES),
53 0 : fUseDefaultTPCPriors(kFALSE)
54 0 : {
55 : //
56 : // default constructor
57 : //
58 0 : for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
59 0 : AliLog::SetClassDebugLevel("AliPIDCombined",10);
60 0 : }
61 :
62 : //-------------------------------------------------------------------------------------------------
63 : AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
64 0 : TNamed(name,title),
65 0 : fDetectorMask(0),
66 0 : fRejectMismatchMask(0x7F),
67 0 : fEnablePriors(kFALSE),
68 0 : fSelectedSpecies(AliPID::kSPECIES),
69 0 : fUseDefaultTPCPriors(kFALSE)
70 0 : {
71 : //
72 : // default constructor with name and title
73 : //
74 0 : for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
75 0 : AliLog::SetClassDebugLevel("AliPIDCombined",10);
76 :
77 0 : }
78 :
79 : //-------------------------------------------------------------------------------------------------
80 0 : AliPIDCombined::~AliPIDCombined() {
81 :
82 0 : for(Int_t i=0;i < AliPID::kSPECIESC;i++){
83 0 : if(fPriorsDistributions[i])
84 0 : delete fPriorsDistributions[i];
85 : }
86 0 : }
87 :
88 : //-------------------------------------------------------------------------------------------------
89 : void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
90 0 : if ( (type < 0) || ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
91 0 : AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
92 0 : return;
93 : }
94 0 : if(prior) {
95 : Int_t i = (Int_t)type;
96 0 : if (fPriorsDistributions[i]) {
97 0 : delete fPriorsDistributions[i];
98 : }
99 0 : fPriorsDistributions[i]=new TH1F(*prior);
100 0 : SetEnablePriors(kTRUE);
101 0 : }
102 0 : }
103 :
104 : //-------------------------------------------------------------------------------------------------
105 : UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities,Double_t *priorsOwn) const {
106 : //
107 : // (1) Get raw probabilities of requested detectors and combine
108 : // (2) priors passed as argument
109 : // (3) Compute Bayes probabilities
110 : //
111 :
112 :
113 : // (1) Get raw probabilities of selected detectors and combine
114 0 : UInt_t usedMask=0; // detectors actually used for track
115 0 : fTOFmismProb = 0; // reset TOF mismatch weights
116 :
117 : AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
118 0 : Double_t p[AliPID::kSPECIESC]; // combined probabilities of selected detectors
119 0 : Double_t pMismTOF[AliPID::kSPECIESC]; // combined TOF mismatch probabilities using selected detectors
120 0 : for (Int_t i=0;i<fSelectedSpecies;i++){ p[i]=1.;pMismTOF[i]=1.;} // no decision
121 0 : for (Int_t ibit = 0; ibit < 7 ; ibit++) {
122 0 : AliPIDResponse::EDetCode detBit = (AliPIDResponse::EDetCode)(1<<ibit);
123 0 : if (fDetectorMask & detBit) { // getting probabilities for requested detectors only
124 0 : Double_t detProb[AliPID::kSPECIESC];
125 0 : status = response->ComputePIDProbability(detBit,track,fSelectedSpecies,detProb);
126 0 : if (status == AliPIDResponse::kDetPidOk) {
127 0 : if (fRejectMismatchMask & detBit) { // mismatch check (currently just for TOF)
128 0 : if (detBit == AliPIDResponse::kDetTOF) {
129 0 : fTOFmismProb = response->GetTOFMismatchProbability(); // mismatch weights computed with TOF probs (no arguments)
130 : //Float_t probMis = response->GetTOFMismatchProbability(track); // mismatch compatibility TPC-TOF cut
131 0 : SetCombinedStatus(status,&usedMask,detBit,detProb,fTOFmismProb);
132 0 : } else {
133 0 : SetCombinedStatus(status,&usedMask,detBit);
134 : }
135 : } else {
136 0 : SetCombinedStatus(status,&usedMask,detBit);
137 : }
138 0 : for (Int_t i=0;i<fSelectedSpecies;i++){
139 0 : p[i] *= detProb[i];
140 0 : if(detBit == AliPIDResponse::kDetTOF){
141 0 : Double_t pt = track->Pt();
142 0 : Double_t mismPropagationFactor[] = {1.,1.,1.,1. + TMath::Exp(1. - 1.12*pt),
143 0 : 1. + 1./(4.71114 - 5.72372*pt + 2.94715*pt*pt),1.,1.,1.,1.,1.}; // correction for kaons and protons which has to be alligned with the one in AliPIDResponse
144 0 : pMismTOF[i] *= fTOFmismProb*mismPropagationFactor[i];
145 0 : }
146 0 : else pMismTOF[i] *= detProb[i];
147 : }
148 0 : }
149 0 : }
150 : }
151 : // if no detectors available there is no point to go further
152 0 : if (usedMask == 0) return usedMask;
153 :
154 : // (2) Get priors and propagate depending on detectors used
155 0 : Double_t priors[AliPID::kSPECIESC];
156 0 : memset(priors,0,fSelectedSpecies*sizeof(Double_t));
157 :
158 0 : if(priorsOwn){
159 0 : for(Int_t ipr=0;ipr < fSelectedSpecies;ipr++)
160 0 : priors[ipr] = priorsOwn[ipr];
161 0 : }
162 :
163 0 : else if (fEnablePriors){
164 0 : Bool_t isPPB = (response->GetBeamType() == AliPIDResponse::kPPB);
165 0 : GetPriors(track,priors,response->GetCurrentCentrality(),isPPB);
166 :
167 : // We apply the propagation factors of the more external detector
168 : //
169 : // TPC+HMPID --> apply HMPID propagation factors (TRD and TOF may be present)
170 : // TPC+EMCAL --> apply EMCAL propagation factors (TRD and TOF may be present)
171 : // TPC+TOF --> apply TOF propagation factors (TRD may be present, HMPID and EMCAL not (if requested))
172 : // TPC+TRD --> apply TRD propagation factors (TOF, HMPID and EMCAL not present (if requested) )
173 : //
174 0 : if(fUseDefaultTPCPriors) {
175 :
176 0 : Double_t pt=TMath::Abs(track->Pt());
177 0 : if ( ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL) || ( (usedMask & AliPIDResponse::kDetHMPID)==AliPIDResponse::kDetHMPID) ) {
178 : // we assume EMCAL and HMPID cannot be simultaneously present
179 0 : if ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL ) {
180 : // EMCal case (for the moment only in combination with TPC)
181 : // propagation factors determined from LHC11d MC (LHC12a15f)
182 : // v2 clusterizer, dEta < 0.015, dPhi < 0.03, NonLinearityFunction = 6
183 :
184 : Float_t electronEMCALfactor=0.1;
185 : Float_t kaonEMCALfactor = 0.1;
186 : Float_t protonEMCALfactor = 0.1;
187 :
188 0 : if(track->Charge() > 0){
189 : // positiv charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
190 0 : if(pt > 0.75 && pt < 20.0){
191 0 : 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)) );
192 0 : 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)) );
193 0 : 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)) );
194 :
195 0 : }
196 : }
197 :
198 0 : if(track->Charge() < 0){
199 : // negative charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
200 0 : if(pt > 0.75 && pt < 20.0){
201 0 : 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)) );
202 0 : 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)) );
203 0 : 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)) );
204 0 : }
205 : }
206 0 : priors[0] *= electronEMCALfactor;
207 0 : priors[3] *= kaonEMCALfactor;
208 0 : priors[4] *= protonEMCALfactor;
209 0 : } // end of EMCAL case
210 :
211 0 : else if ( (usedMask & AliPIDResponse::kDetHMPID)==AliPIDResponse::kDetHMPID ) { // HMPID case
212 :
213 : Float_t kaonHMPIDfactor = 0.;
214 : Float_t protonHMPIDfactor = 0.;
215 0 : if(pt>1. && pt<6.) kaonHMPIDfactor = (-0.0729337 + pt*0.0999531 - pt*pt*0.0371803 + pt*pt*pt*0.00706436 - pt*pt*pt*pt*0.000643619 + pt*pt*pt*pt*pt*2.21853e-05)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06);
216 0 : if(pt>1.4 && pt<6.) protonHMPIDfactor = (-0.0444188 + pt*0.0681506 - pt*pt*0.0231819 + pt*pt*pt*0.00400771 - pt*pt*pt*pt*0.000339315 + pt*pt*pt*pt*pt*1.12616e-05)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06);
217 0 : if(pt>6. && pt<8.) kaonHMPIDfactor = (-0.0729337 + pt*0.0999531 - pt*pt*0.0371803 + pt*pt*pt*0.00706436 - pt*pt*pt*pt*0.000643619 + pt*pt*pt*pt*pt*2.21853e-05)/0.0530456;
218 0 : if(pt>8.) kaonHMPIDfactor = 0.0550432/0.0530456;
219 0 : if(pt>6. && pt<8.5) protonHMPIDfactor = (-0.0444188 + pt*0.0681506 - pt*pt*0.0231819 + pt*pt*pt*0.00400771 - pt*pt*pt*pt*0.000339315 + pt*pt*pt*pt*pt*1.12616e-05)/0.0530456;
220 0 : if(pt>8.5) protonHMPIDfactor = 0.0530071/0.0530456;
221 :
222 0 : if(track->Charge() < 0){
223 0 : if(pt>0.4 && pt<6.) protonHMPIDfactor = (-0.0351485 + pt*0.0473821 - pt*pt*0.0147947 + pt*pt*pt*0.00254811- pt*pt*pt*pt*0.000224724 + pt*pt*pt*pt*pt*7.9303e-06)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06);
224 0 : if(pt>6. && pt<8.5) protonHMPIDfactor = (-0.0351485 + pt*0.0473821 - pt*pt*0.0147947 + pt*pt*pt*0.00254811- pt*pt*pt*pt*0.000224724 + pt*pt*pt*pt*pt*7.9303e-06)/0.0530456;
225 0 : if(pt>8.5) protonHMPIDfactor = 0.0457756/0.0530456;
226 : }
227 :
228 0 : priors[3] *= kaonHMPIDfactor;
229 0 : priors[4] *= protonHMPIDfactor;
230 :
231 0 : }
232 :
233 : } // end of outer cases: EMCAL/HMPID
234 0 : else if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){
235 : Float_t kaonTOFfactor = 0.1;
236 0 : if(pt > 0.29){
237 0 : kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/4.96502e-01)*TMath::Power(pt,-1.50705);
238 0 : }
239 :
240 : Float_t protonTOFfactor = 0.1;
241 0 : if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
242 :
243 0 : if(track->Charge() < 0){ // for negative tracks
244 0 : kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
245 0 : protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
246 0 : }
247 : // TODO: we may need an electron factor for TOF as well, especially if TRD is present!
248 0 : priors[3] *= kaonTOFfactor;
249 0 : priors[4] *= protonTOFfactor;
250 0 : } // end of TOF case
251 0 : else if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
252 : Float_t electronTRDfactor=0.1;
253 : Float_t kaonTRDfactor = 0.1;
254 : Float_t protonTRDfactor = 0.1;
255 :
256 0 : if(track->Charge() > 0){
257 : // positiv charge
258 0 : if(pt > 0.25) electronTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
259 0 : if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
260 0 : if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
261 : }
262 :
263 0 : if(track->Charge() < 0){
264 : // negative charge
265 0 : if(pt > 0.25) electronTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
266 0 : if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
267 0 : if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
268 : }
269 : // what about electrons
270 0 : priors[0] *= electronTRDfactor;
271 0 : priors[3] *= kaonTRDfactor;
272 0 : priors[4] *= protonTRDfactor;
273 0 : } // end of TRD case
274 0 : } // end of fUseDefaultTPCPriors
275 0 : } // end of use priors
276 0 : else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
277 :
278 : // Compute Bayes probabilities
279 0 : ComputeBayesProbabilities(bayesProbabilities,p,priors,pMismTOF);
280 :
281 : // compute TOF probability contribution from mismatch
282 0 : fTOFmismProb = 0;
283 0 : for (Int_t i=0;i<fSelectedSpecies;i++) fTOFmismProb += pMismTOF[i];
284 :
285 0 : return usedMask;
286 0 : }
287 :
288 :
289 : //-------------------------------------------------------------------------------------------------
290 : void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality,Bool_t isPPB) const {
291 :
292 : //
293 : // get priors from distributions
294 : //
295 :
296 :
297 0 : Double_t pt=TMath::Abs(track->Pt());
298 :
299 0 : if(fUseDefaultTPCPriors){ // use default priors if requested
300 0 : Float_t usedCentr = centrality+5*(centrality>0);
301 0 : if(usedCentr < -0.99) usedCentr = -0.99;
302 0 : else if(usedCentr > 99) usedCentr = 99;
303 0 : if(pt > 9.99) pt = 9.99;
304 0 : else if(pt < 0.01) pt = 0.01;
305 :
306 0 : if(! isPPB)
307 0 : for(Int_t i=0;i<AliPID::kSPECIESC;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
308 : else
309 0 : for(Int_t i=0;i<AliPID::kSPECIESC;i++) p[i] = fDefaultPriorsTPCpPb[i]->Interpolate(usedCentr,pt);
310 :
311 : return;
312 : }
313 :
314 : Double_t sumPriors = 0;
315 0 : for (Int_t i=0;i<fSelectedSpecies;++i){
316 0 : if (fPriorsDistributions[i]){
317 0 : p[i]=fPriorsDistributions[i]->Interpolate(pt);
318 0 : }
319 : else {
320 0 : p[i]=0.;
321 : }
322 0 : sumPriors+=p[i];
323 : }
324 :
325 : // normalizing priors
326 0 : if (sumPriors == 0) return;
327 0 : for (Int_t i=0;i<fSelectedSpecies;++i){
328 0 : p[i]=p[i]/sumPriors;
329 : }
330 0 : return;
331 0 : }
332 : //-------------------------------------------------------------------------------------------------
333 : void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
334 :
335 : //
336 : // get priors from distributions
337 : //
338 0 : Double_t pt=TMath::Abs(track->Pt());
339 :
340 0 : if(fUseDefaultTPCPriors){ // use default priors if requested
341 0 : Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
342 0 : if(usedCentr < -0.99) usedCentr = -0.99;
343 0 : else if(usedCentr > 99) usedCentr = 99;
344 0 : if(pt > 9.99) pt = 9.99;
345 0 : else if(pt < 0.01) pt = 0.01;
346 :
347 0 : for(Int_t i=0;i<AliPID::kSPECIESC;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
348 :
349 : // Extra factor if TOF matching was required
350 0 : if(detUsed & AliPIDResponse::kDetTOF){
351 : Float_t kaonTOFfactor = 0.1;
352 0 : if(pt > 0.29){
353 0 : kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/4.96502e-01)*TMath::Power(pt,-1.50705);
354 0 : }
355 : Float_t protonTOFfactor = 0.1;
356 0 : if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
357 :
358 0 : if(track->Charge() < 0){ // for negative tracks
359 0 : kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
360 0 : protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
361 0 : }
362 :
363 0 : p[3] *= kaonTOFfactor;
364 0 : p[4] *= protonTOFfactor;
365 0 : }
366 :
367 : return;
368 : }
369 :
370 :
371 : Double_t sumPriors = 0;
372 0 : for (Int_t i=0;i<fSelectedSpecies;++i){
373 0 : if (fPriorsDistributions[i]){
374 0 : p[i]=fPriorsDistributions[i]->Interpolate(pt);
375 0 : }
376 : else {
377 0 : p[i]=0.;
378 : }
379 0 : sumPriors+=p[i];
380 : }
381 :
382 : // normalizing priors
383 0 : if (sumPriors == 0) return;
384 0 : for (Int_t i=0;i<fSelectedSpecies;++i){
385 0 : p[i]=p[i]/sumPriors;
386 : }
387 0 : return;
388 0 : }
389 : //-------------------------------------------------------------------------------------------------
390 : void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior, Double_t* probDensityMism) const {
391 :
392 :
393 : //
394 : // calculate Bayesian probabilities
395 : //
396 : Double_t sum = 0.;
397 0 : for (Int_t i = 0; i < fSelectedSpecies; i++) {
398 0 : sum += probDensity[i] * prior[i];
399 : }
400 0 : if (sum <= 0) {
401 :
402 0 : AliError("Invalid probability densities or priors");
403 0 : for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = 1./fSelectedSpecies;
404 0 : return;
405 : }
406 0 : for (Int_t i = 0; i < fSelectedSpecies; i++) {
407 0 : probabilities[i] = probDensity[i] * prior[i] / sum;
408 0 : if(probDensityMism) probDensityMism[i] *= prior[i] / sum;
409 : }
410 :
411 :
412 0 : }
413 :
414 :
415 : //----------------------------------------------------------------------------------------
416 : void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit) const {
417 0 : switch (status) {
418 : case AliPIDResponse::kDetNoParams:
419 : case AliPIDResponse::kDetNoSignal:
420 : case AliPIDResponse::kDetMismatch: // for backward compatibilty, we need then to remove kDetMismatch from AliPIDResponse
421 : break;
422 : case AliPIDResponse::kDetPidOk:
423 0 : *maskDetIn = *maskDetIn | bit;
424 0 : break;
425 : }
426 :
427 0 : }
428 :
429 : //----------------------------------------------------------------------------------------
430 : void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* /*p*/, const Float_t /*probMis*/) const {
431 0 : switch (status) {
432 : case AliPIDResponse::kDetNoParams:
433 : case AliPIDResponse::kDetNoSignal:
434 : case AliPIDResponse::kDetMismatch: // for backward compatibility, we need then to remove kDeteMismatch from AliPIDResponse
435 : break;
436 : case AliPIDResponse::kDetPidOk:
437 : // if (probMis > 0.5) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies; // no longer used because mismatch is in the framework now
438 : //else
439 0 : *maskDetIn = *maskDetIn | bit;
440 0 : break;
441 : }
442 :
443 0 : }
444 :
445 :
446 :
447 : //----------------------------------------------------------------------------------------
448 : void AliPIDCombined::SetDefaultTPCPriors(){
449 0 : fEnablePriors=kTRUE;
450 0 : fUseDefaultTPCPriors = kTRUE;
451 :
452 : //Check if priors are already initialized
453 0 : if (fDefaultPriorsTPC[0]) return;
454 :
455 0 : TString oadbfilename("$ALICE_PHYSICS/OADB/COMMON/PID/data/");
456 : // TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
457 : // TString oadbfilename(Form("%s/COMMON/PID/data/",AliAnalysisManager::GetOADBPath()));
458 0 : oadbfilename += "/PIDdefaultPriors.root";
459 0 : TFile * foadb = TFile::Open(oadbfilename.Data());
460 0 : if(!foadb || !foadb->IsOpen()) {
461 0 : delete foadb;
462 0 : AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
463 0 : return;
464 : }
465 :
466 0 : AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
467 0 : if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
468 :
469 0 : const char *namespecies[AliPID::kSPECIESC] = {"El","Mu","Pi","Ka","Pr","De","Tr","He3","He3"};
470 0 : char name[100];
471 :
472 0 : for(Int_t i=0;i < AliPID::kSPECIESC;i++){
473 0 : snprintf(name,99,"hDefault%sPriors",namespecies[i]);
474 0 : fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
475 0 : if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
476 0 : snprintf(name,99,"hDefault%sPriorsPPb",namespecies[i]);
477 0 : fDefaultPriorsTPCpPb[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
478 0 : if (!fDefaultPriorsTPCpPb[i]) {
479 0 : fDefaultPriorsTPCpPb[i] = fDefaultPriorsTPC[i]; // use PbPb ones
480 0 : }
481 : }
482 :
483 0 : delete foadb;
484 0 : }
|
