9dc68f2355956aefaee50eef1e802bef572f6f94
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAODPidHF.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  * *************************************************************************/
16
17 //***********************************************************
18 // Class AliAODPidHF
19 // class for PID with AliAODRecoDecayHF
20 // Authors: D. Caffarri caffarri@pd.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de
21 //***********************************************************
22 #include "AliAODPidHF.h"
23 #include "AliAODPid.h"
24 #include "AliAODTrack.h"
25 #include "AliPID.h"
26 #include "AliTPCPIDResponse.h"
27 #include "AliITSPIDResponse.h"
28 #include "AliTOFPIDResponse.h"
29 #include "AliAODpidUtil.h"
30 #include "AliESDtrack.h"
31
32
33 ClassImp(AliAODPidHF)
34
35 //------------------------------
36 AliAODPidHF::AliAODPidHF():
37   AliAODPid(),
38   fnNSigma(5),
39   fnSigma(),
40   fTOFSigma(160.),
41   fnPriors(5),
42   fPriors(),
43   fnPLimit(2),
44   fPLimit(),
45   fAsym(kFALSE)
46 {
47  //
48  // Default constructor
49  //
50  fPLimit=new Double_t[fnPLimit];
51  fnSigma=new Double_t[fnNSigma];
52  fPriors=new Double_t[fnPriors];
53
54  for(Int_t i=0;i<fnNSigma;i++){
55   fnSigma[i]=0.;
56  }
57  for(Int_t i=0;i<fnPriors;i++){
58   fPriors[i]=0.;
59  }
60  for(Int_t i=0;i<fnPLimit;i++){
61   fPLimit[i]=0.;
62  }
63
64 }
65 //----------------------
66 AliAODPidHF::~AliAODPidHF()
67 {
68       // destructor
69  //   if(fPLimit) delete fPLimit;
70  //   if(fnSigma)  delete fnSigma;
71  //   if(fPriors)  delete fPriors;
72 }
73 //------------------------
74 AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
75   AliAODPid(pid),
76   fnNSigma(pid.fnNSigma),
77   fnSigma(pid.fnSigma),
78   fTOFSigma(pid.fTOFSigma),
79   fnPriors(pid.fnPriors),
80   fPriors(pid.fPriors),
81   fnPLimit(pid.fnPLimit),
82   fPLimit(pid.fPLimit),
83   fAsym(pid.fAsym)
84   {
85   
86   for(Int_t i=0;i<5;i++){
87     fPriors[i]=pid.fPriors[i];
88   }
89   
90   }
91
92 //----------------------
93 Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
94 // raw PID for single detectors, returns the particle type with smaller sigma
95    Int_t specie=-1;
96    if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
97    if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
98    if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
99
100   return specie;
101
102 }
103 //---------------------------
104 Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
105 // checks if the track can be a kaon, raw PID applied for single detectors
106  Int_t specie=0;
107
108  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
109  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
110  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
111
112  if(specie==3) return kTRUE;
113  return kFALSE;
114 }
115 //---------------------------
116 Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
117 // checks if the track can be a pion, raw PID applied for single detectors
118
119  Int_t specie=0;
120
121  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
122  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
123  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
124
125  if(specie==2) return kTRUE;
126  return kFALSE;
127 }
128 //---------------------------
129 Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
130 // checks if the track can be a proton raw PID applied for single detectors
131
132  Int_t specie=0;
133  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
134  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4); 
135  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
136
137  if(specie==4) return kTRUE;
138
139  return kFALSE;
140 }
141 //--------------------------
142 Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
143 // checks if the track can be an electron raw PID applied for single detectors
144
145  Int_t specie=-1;
146  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
147  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
148  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
149
150  if(specie==0) return kTRUE;
151
152  return kFALSE;
153 }
154 //--------------------------
155 Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
156 // n-sigma cut, TPC PID
157
158   if(!CheckStatus(track,"TPC")) return -1;
159   AliAODPid *pidObj = track->GetDetPid();
160   
161   Double_t dedx=pidObj->GetTPCsignal();
162   Double_t mom = pidObj->GetTPCmomentum();
163   AliTPCPIDResponse tpcResponse;
164   Int_t pid=-1;
165   if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
166    Double_t nsigmaMax=fnSigma[0];
167    for(Int_t ipart=0;ipart<5;ipart++){
168     AliPID::EParticleType type=AliPID::EParticleType(ipart);
169     Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
170     if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
171      pid=ipart;
172      nsigmaMax=nsigma;
173     }
174    }
175   }else{ // asks only for one particle specie
176    AliPID::EParticleType type=AliPID::EParticleType(specie);
177     Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
178    if (nsigma>fnSigma[0]) {
179     pid=-1; 
180    }else{
181     pid=specie;
182    }
183   }
184
185  return pid;
186
187 }
188 //----------------------------
189 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
190 // truncated mean, ITS PID
191
192   if(!CheckStatus(track,"ITS")) return -1;
193
194   Double_t mom=track->P();
195   AliAODPid *pidObj = track->GetDetPid();
196
197   Double_t dedx=pidObj->GetITSsignal();
198   AliITSPIDResponse itsResponse;
199   Int_t pid=-1;
200   if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
201    Double_t nsigmaMax=fnSigma[4];
202    for(Int_t ipart=0;ipart<5;ipart++){
203     AliPID::EParticleType type=AliPID::EParticleType(ipart);
204     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
205     if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
206      pid=ipart;
207      nsigmaMax=nsigma;
208     }
209    }
210   }else{ // asks only for one particle specie
211    AliPID::EParticleType type=AliPID::EParticleType(specie);
212     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
213    if (nsigma>fnSigma[4]) {
214     pid=-1; 
215    }else{
216     pid=specie;
217    }
218   }
219  return pid; 
220 }
221 //----------------------------
222 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
223 // n-sigma cut, TOF PID
224
225  if(!CheckStatus(track,"TOF")) return -1;
226
227  Double_t time[AliPID::kSPECIESN];
228  AliAODPid *pidObj = track->GetDetPid();
229  pidObj->GetIntegratedTimes(time);
230  Double_t sigTOF=pidObj->GetTOFsignal();
231 // AliTOFPIDResponse tofResponse;
232  Int_t pid=-1;
233
234   if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
235    Double_t nsigmaMax=fTOFSigma*fnSigma[3];
236    for(Int_t ipart=0;ipart<5;ipart++){
237     //AliPID::EParticleType type=AliPID::EParticleType(ipart);
238     //Double_t nsigma = tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type));
239     Double_t nsigma=TMath::Abs(sigTOF-time[ipart]);
240     if((nsigma<nsigmaMax) && (nsigma<fnSigma[3]*fTOFSigma)) {
241      pid=ipart;
242      nsigmaMax=nsigma;
243     }
244    }
245   }else{ // asks only for one particle specie
246    //AliPID::EParticleType type=AliPID::EParticleType(specie);
247    //Double_t nsigma = TMath::Abs(tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type)));
248     Double_t nsigma=TMath::Abs(sigTOF-time[specie]);
249    if (nsigma>fnSigma[3]*fTOFSigma) {
250     pid=-1; 
251    }else{
252     pid=specie;
253    }
254   }
255  return pid; 
256
257 }
258 //------------------------------
259 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
260 // combined PID stored inside the AOD track
261
262  const Double_t *pid=track->PID();
263  Float_t max=0.;
264  Int_t k=-1;
265  for (Int_t i=0; i<10; i++) {
266    if (pid[i]>max) {k=i; max=pid[i];}  
267  }
268
269  if(k==2) type[0]=kTRUE;
270  if(k==3) type[1]=kTRUE;
271  if(k==4) type[2]=kTRUE;
272
273  return;
274 }
275 //--------------------
276 void AliAODPidHF::BayesianProbability(AliAODTrack *track,TString detectors,Double_t *pid) const{
277 // bayesian PID for single detectors or combined
278
279   if(detectors.Contains("ITS")) {BayesianProbabilityITS(track,pid);return;}
280   if(detectors.Contains("TPC")) {BayesianProbabilityTPC(track,pid);return;}
281   if(detectors.Contains("TOF")) {BayesianProbabilityTOF(track,pid);return;}
282
283   if(detectors.Contains("All")) {
284     Double_t probITS[5]={0.,0.,0.,0.,0.};
285     Double_t probTPC[5]={0.,0.,0.,0.,0.};
286     Double_t probTOF[5]={0.,0.,0.,0.,0.};
287     BayesianProbabilityITS(track,probITS);
288     BayesianProbabilityTPC(track,probTPC);
289     BayesianProbabilityTOF(track,probTOF);
290     Double_t probTot[5]={0.,0.,0.,0.,0.};
291     for(Int_t i=0;i<5;i++){
292      probTot[i]=probITS[i]*probTPC[i]*probTOF[i];
293     }
294     for(Int_t i2=0;i2<5;i2++){
295      pid[i2]=probTot[i2]*fPriors[i2]/(probTot[0]*fPriors[0]+probTot[1]*fPriors[1]+probTot[2]*fPriors[2]+probTot[3]*fPriors[3]+probTot[4]*fPriors[4]);
296     }
297    }
298
299  return;
300
301 }
302 //------------------------------------
303 void AliAODPidHF::BayesianProbabilityITS(AliAODTrack *track,Double_t *prob) const{
304
305 // bayesian PID for ITS
306  AliAODpidUtil pid;
307  Double_t itspid[AliPID::kSPECIES];
308  pid.MakeITSPID(track,itspid);
309  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
310   prob[ind]=itspid[ind]*fPriors[ind]/(itspid[0]*fPriors[0]+itspid[1]*fPriors[1]+itspid[2]*fPriors[2]+itspid[3]*fPriors[3]+itspid[4]*fPriors[4]);
311  }
312  return;
313
314 }
315 //------------------------------------
316 void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack *track,Double_t *prob) const{
317 // bayesian PID for TPC
318
319  AliAODpidUtil pid;
320  Double_t tpcpid[AliPID::kSPECIES];
321  pid.MakeTPCPID(track,tpcpid);
322  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
323   if(tpcpid[ind]>0.) {
324    prob[ind]=tpcpid[ind]*fPriors[ind]/(tpcpid[0]*fPriors[0]+tpcpid[1]*fPriors[1]+tpcpid[2]*fPriors[2]+tpcpid[3]*fPriors[3]+tpcpid[4]*fPriors[4]);
325   }else{
326    prob[ind]=0.;
327   }
328  }
329  return;
330
331 }
332 //------------------------------------
333 void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack *track,Double_t *prob) const{
334 // bayesian PID for TOF
335
336  AliAODpidUtil pid;
337  Double_t tofpid[AliPID::kSPECIES];
338  pid.MakeTOFPID(track,tofpid);
339  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
340   prob[ind]=tofpid[ind]*fPriors[ind]/(tofpid[0]*fPriors[0]+tofpid[1]*fPriors[1]+tofpid[2]*fPriors[2]+tofpid[3]*fPriors[3]+tofpid[4]*fPriors[4]);
341  }
342  return;
343
344 }
345 //---------------------------------
346 void AliAODPidHF::BayesianProbabilityTRD(AliAODTrack *track,Double_t *prob) const{
347 // bayesian PID for TRD
348
349  AliAODpidUtil pid;
350  Double_t trdpid[AliPID::kSPECIES];
351  pid.MakeTRDPID(track,trdpid);
352  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
353   if(trdpid[ind]>0.) {
354    prob[ind]=trdpid[ind]*fPriors[ind]/(trdpid[0]*fPriors[0]+trdpid[1]*fPriors[1]+trdpid[2]*fPriors[2]+trdpid[3]*fPriors[3]+trdpid[4]*fPriors[4]);
355   }else{
356    prob[ind]=0.;
357   }
358  }
359   return;
360
361  }
362 //--------------------------------
363 Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
364
365 // Quality cuts on the tracks, detector by detector
366
367  if(detectors.Contains("ITS")){
368   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
369   UChar_t clumap=track->GetITSClusterMap();
370   Int_t nPointsForPid=0;
371   for(Int_t i=2; i<6; i++){
372    if(clumap&(1<<i)) ++nPointsForPid;
373   }
374   if(nPointsForPid<3) return kFALSE;// track not to be used for PID purposes
375  }
376
377  if(detectors.Contains("TPC")){
378    if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
379    UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
380    if (nTPCClus<70) return kFALSE;
381  }
382
383  if(detectors.Contains("TOF")){
384    if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return kFALSE;
385    if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return kFALSE;
386    if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return kFALSE;
387  }
388
389
390  if(detectors.Contains("TRD")){
391   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return kFALSE;
392   AliAODPid *pidObj = track->GetDetPid();
393   Float_t *mom=pidObj->GetTRDmomentum();
394   Int_t ntracklets=0;
395   for(Int_t iPl=0;iPl<6;iPl++){
396    if(mom[iPl]>0.) ntracklets++;
397   }
398    if(ntracklets<4) return kFALSE;
399  }
400
401  return kTRUE;
402 }
403 //--------------------------------------------
404 Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
405 // TPC nsigma cut PID, different sigmas in different p bins
406
407   if(!CheckStatus(track,"TPC")) return kFALSE;
408   AliAODPid *pidObj = track->GetDetPid();
409   Double_t mom = pidObj->GetTPCmomentum();
410   Double_t dedx=pidObj->GetTPCsignal();
411   
412   AliTPCPIDResponse tpcResponse;
413   AliPID::EParticleType type=AliPID::EParticleType(specie);
414   Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type)); 
415
416   if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
417   if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
418   if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
419
420  return kFALSE;
421 }
422 //------------------
423 Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track,Int_t mode,Int_t specie,Bool_t compat){
424 // combination of the PID info coming from TPC and TOF
425  if(mode==1){
426   //TOF || TPC (a la' Andrea R.)
427  // convention (temporary): 
428  // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
429  // the method returns the sum of the response of the 2 detectors
430   if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;
431
432   Int_t TPCinfo=0;
433   if(CheckStatus(track,"TPC")) {
434    if(fAsym) {
435     if(TPCRawAsym(track,specie)) {
436       TPCinfo=1;
437      }else{
438       TPCinfo=-1;
439      }
440    }else{
441     if(specie==2 && IsPionRaw(track,"TPC")) {
442      TPCinfo=1;
443     }else{
444      TPCinfo=-1;
445     }
446     if(specie==3 && IsKaonRaw(track,"TPC")) {
447      TPCinfo=1;
448     }else{
449      TPCinfo=-1;
450     }
451     if(specie==4 && IsProtonRaw(track,"TPC")) {
452      TPCinfo=1;
453     }else{
454      TPCinfo=-1;
455     }
456
457    }
458
459
460   if(compat && TPCinfo<0){
461    Double_t sig0tmp=fnSigma[0];
462    SetSigma(0,3.);
463    if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=0;
464    if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=0;
465    if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=0;
466    SetSigma(0,sig0tmp);
467   }
468
469  }
470
471  if(!CheckStatus(track,"TOF")) return TPCinfo;
472
473  Int_t TOFinfo=-1;
474  
475   if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
476   if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
477   if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;
478
479  if(compat && TOFinfo>0){
480   Double_t ptrack=track->P();
481   if(ptrack>1.5) TOFinfo=0;
482  }
483
484  return TPCinfo+TOFinfo;
485  }
486  if(mode==2){
487   //TPC & TOF (a la' Yifei)
488  // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
489   Int_t TPCinfo=1; 
490   if(CheckStatus(track,"TPC")) {
491    if(fAsym){
492     if(!TPCRawAsym(track,specie)) TPCinfo=-1;
493    }else{
494     if(specie==2 && !IsPionRaw(track,"TPC")) TPCinfo=-1;
495     if(specie==3 && !IsKaonRaw(track,"TPC")) TPCinfo=-1;
496     if(specie==4 && !IsProtonRaw(track,"TPC")) TPCinfo=-1;
497    }
498   }
499
500   Int_t TOFinfo=1;
501   if(!CheckStatus(track,"TOF")) return TPCinfo;
502
503    if(specie==2 && !IsPionRaw(track,"TOF")) TOFinfo=-1;
504    if(specie==3 && !IsKaonRaw(track,"TOF")) TOFinfo=-1;
505    if(specie==4 && !IsProtonRaw(track,"TOF")) TOFinfo=-1;
506
507    if(TOFinfo==1 && TPCinfo==1) return 1;
508    return -1;
509
510  }
511
512  if(mode==3){
513  //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
514  // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
515   
516   if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;
517
518   Double_t ptrack=track->P();
519   
520
521   Int_t TPCinfo=-1;
522
523    if(ptrack<fPLimit[0]) {  
524     if(!CheckStatus(track,"TPC")) return 0;
525     if(fAsym) {
526      if(TPCRawAsym(track,specie)) TPCinfo=1;
527     }else{
528      if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=1;
529      if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=1;
530      if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=1;
531     } 
532     return TPCinfo;
533    }
534
535    Int_t TOFinfo=-1;
536    if(ptrack>=fPLimit[0]){
537     if(!CheckStatus(track,"TOF")) return 0;
538     if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
539     if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
540     if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;
541     return TOFinfo;
542    }
543
544  }
545
546  return -1;
547 }
548