]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAODPidHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Speed-up code for D meson PID
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAODPidHF.cxx
CommitLineData
7ce8802c 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
27de2dfb 17/* $Id$ */
18
7ce8802c 19//***********************************************************
20// Class AliAODPidHF
21// class for PID with AliAODRecoDecayHF
b257e565 22// 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 P. Antonioli pietro.antonioli@bo.infn.it
7ce8802c 23//***********************************************************
64e39959 24#include <TCanvas.h>
25
7ce8802c 26#include "AliAODPidHF.h"
27#include "AliAODPid.h"
7ce8802c 28#include "AliPID.h"
b79bfc3e 29#include "AliPIDResponse.h"
7ce8802c 30#include "AliAODpidUtil.h"
f9e3cf8e 31#include "AliESDtrack.h"
7ce8802c 32
33
34ClassImp(AliAODPidHF)
35
36//------------------------------
37AliAODPidHF::AliAODPidHF():
228f3aba 38 AliAODPid(),
39 fnNSigma(5),
40 fnSigma(),
41 fTOFSigma(160.),
42 fnPriors(5),
43 fPriors(),
44 fnPLimit(2),
45 fPLimit(),
cdac45d3 46 fAsym(kFALSE),
47 fTPC(kFALSE),
48 fTOF(kFALSE),
49 fITS(kFALSE),
50 fTRD(kFALSE),
79f3c225 51 fMatch(0),
913d1957 52 fCompat(kFALSE),
ad42e35b 53 fPCompatTOF(1.5),
02baac36 54 fnNSigmaCompat(2),
ad42e35b 55 fnSigmaCompat(),
12141139 56 fMC(kFALSE),
17c4a8d5 57 fOnePad(kFALSE),
5821ea5b 58 fMCLowEn2011(kFALSE),
67a82ef2 59 fppLowEn2011(kFALSE),
439d61b7 60 fPbPb(kFALSE),
b79bfc3e 61 fTOFdecide(kFALSE),
2c9eefd4 62 fOldPid(kTRUE),
63 fPtThresholdTPC(999999.),
11690a06 64 fPidResponse(0),
70e398d3 65 fPidCombined(new AliPIDCombined()),
66 fTPCResponse(new AliTPCPIDResponse())
7ce8802c 67{
68 //
69 // Default constructor
70 //
228f3aba 71 fPLimit=new Double_t[fnPLimit];
72 fnSigma=new Double_t[fnNSigma];
73 fPriors=new Double_t[fnPriors];
02baac36 74 fnSigmaCompat=new Double_t[fnNSigmaCompat];
228f3aba 75
76 for(Int_t i=0;i<fnNSigma;i++){
77 fnSigma[i]=0.;
78 }
79 for(Int_t i=0;i<fnPriors;i++){
80 fPriors[i]=0.;
81 }
82 for(Int_t i=0;i<fnPLimit;i++){
83 fPLimit[i]=0.;
84 }
02baac36 85 for(Int_t i=0;i<fnNSigmaCompat;i++){
ad42e35b 86 fnSigmaCompat[i]=3.;
87 }
7ce8802c 88
89}
90//----------------------
91AliAODPidHF::~AliAODPidHF()
92{
93 // destructor
228f3aba 94 // if(fPLimit) delete fPLimit;
95 // if(fnSigma) delete fnSigma;
96 // if(fPriors) delete fPriors;
70e398d3 97 delete fTPCResponse;
7ce8802c 98}
99//------------------------
100AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
101 AliAODPid(pid),
228f3aba 102 fnNSigma(pid.fnNSigma),
375df9ce 103 fnSigma(0),
228f3aba 104 fTOFSigma(pid.fTOFSigma),
105 fnPriors(pid.fnPriors),
375df9ce 106 fPriors(0),
228f3aba 107 fnPLimit(pid.fnPLimit),
375df9ce 108 fPLimit(0),
cdac45d3 109 fAsym(pid.fAsym),
110 fTPC(pid.fTPC),
111 fTOF(pid.fTOF),
112 fITS(pid.fITS),
113 fTRD(pid.fTRD),
114 fMatch(pid.fMatch),
913d1957 115 fCompat(pid.fCompat),
ad42e35b 116 fPCompatTOF(pid.fPCompatTOF),
02baac36 117 fnNSigmaCompat(pid.fnNSigmaCompat),
ad42e35b 118 fnSigmaCompat(pid.fnSigmaCompat),
12141139 119 fMC(pid.fMC),
17c4a8d5 120 fOnePad(pid.fOnePad),
5821ea5b 121 fMCLowEn2011(pid.fMCLowEn2011),
67a82ef2 122 fppLowEn2011(pid.fppLowEn2011),
439d61b7 123 fPbPb(pid.fPbPb),
b79bfc3e 124 fTOFdecide(pid.fTOFdecide),
125 fOldPid(pid.fOldPid),
2c9eefd4 126 fPtThresholdTPC(pid.fPtThresholdTPC),
11690a06 127 fPidResponse(pid.fPidResponse),
70e398d3 128 fPidCombined(pid.fPidCombined),
129 fTPCResponse(pid.fTPCResponse)
375df9ce 130{
7ce8802c 131
375df9ce 132 fnSigma = new Double_t[fnNSigma];
133 for(Int_t i=0;i<fnNSigma;i++){
134 fnSigma[i]=pid.fnSigma[i];
135 }
136 fPriors = new Double_t[fnPriors];
137 for(Int_t i=0;i<fnPriors;i++){
7ce8802c 138 fPriors[i]=pid.fPriors[i];
139 }
375df9ce 140 fPLimit = new Double_t[fnPLimit];
141 for(Int_t i=0;i<fnPLimit;i++){
142 fPLimit[i]=pid.fPLimit[i];
7ce8802c 143 }
144
375df9ce 145
146 //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
147 //fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
148
149}
7ce8802c 150//----------------------
228f3aba 151Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
152// raw PID for single detectors, returns the particle type with smaller sigma
f9e3cf8e 153 Int_t specie=-1;
154 if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
155 if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
156 if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
7ce8802c 157
f9e3cf8e 158 return specie;
7ce8802c 159
160}
161//---------------------------
228f3aba 162Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
163// checks if the track can be a kaon, raw PID applied for single detectors
7ce8802c 164 Int_t specie=0;
165
166 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
167 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
168 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
169
170 if(specie==3) return kTRUE;
171 return kFALSE;
172}
173//---------------------------
228f3aba 174Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
175// checks if the track can be a pion, raw PID applied for single detectors
7ce8802c 176
177 Int_t specie=0;
178
179 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
180 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
181 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
182
183 if(specie==2) return kTRUE;
184 return kFALSE;
185}
186//---------------------------
228f3aba 187Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
188// checks if the track can be a proton raw PID applied for single detectors
7ce8802c 189
190 Int_t specie=0;
191 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
192 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
193 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
194
195 if(specie==4) return kTRUE;
196
197 return kFALSE;
198}
f9e3cf8e 199//--------------------------
228f3aba 200Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
201// checks if the track can be an electron raw PID applied for single detectors
f9e3cf8e 202
203 Int_t specie=-1;
204 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
205 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
206 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
207
208 if(specie==0) return kTRUE;
209
210 return kFALSE;
211}
212//--------------------------
228f3aba 213Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
214// n-sigma cut, TPC PID
7ce8802c 215
70e398d3 216 if(!CheckTPCPIDStatus(track)) return 0;
b79bfc3e 217 Int_t pid=-1;
218 if(fOldPid){
7ce8802c 219 AliAODPid *pidObj = track->GetDetPid();
220
221 Double_t dedx=pidObj->GetTPCsignal();
222 Double_t mom = pidObj->GetTPCmomentum();
2c9eefd4 223 if(mom>fPtThresholdTPC) return 0;
0aba57ed 224 UShort_t nTPCClus=pidObj->GetTPCsignalN();
225 if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
12141139 226
f9e3cf8e 227 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
228f3aba 228 Double_t nsigmaMax=fnSigma[0];
f9e3cf8e 229 for(Int_t ipart=0;ipart<5;ipart++){
7ce8802c 230 AliPID::EParticleType type=AliPID::EParticleType(ipart);
70e398d3 231 Double_t nsigma = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
228f3aba 232 if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
7ce8802c 233 pid=ipart;
234 nsigmaMax=nsigma;
235 }
236 }
237 }else{ // asks only for one particle specie
238 AliPID::EParticleType type=AliPID::EParticleType(specie);
70e398d3 239 Double_t nsigma = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
228f3aba 240 if (nsigma>fnSigma[0]) {
f9e3cf8e 241 pid=-1;
7ce8802c 242 }else{
79f3c225 243 pid=specie;
7ce8802c 244 }
245 }
b79bfc3e 246 }else{ //old pid
247 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
248 Double_t nsigmaMax=fnSigma[0];
249 for(Int_t ipart=0;ipart<5;ipart++){
250 AliPID::EParticleType type=AliPID::EParticleType(ipart);
251 Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
252 if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
253 pid=ipart;
254 nsigmaMax=nsigma;
255 }
256 }
257 }else{ // asks only for one particle specie
258 AliPID::EParticleType type=AliPID::EParticleType(specie);
259 Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
260 if (nsigma>fnSigma[0]) {
261 pid=-1;
262 }else{
263 pid=specie;
264 }
265 }
266
267 } //new pid
7ce8802c 268
269 return pid;
270
271}
272//----------------------------
228f3aba 273Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
274// truncated mean, ITS PID
f9e3cf8e 275
70e398d3 276 if(!CheckITSPIDStatus(track)) return 0;
b79bfc3e 277 Int_t pid=-1;
f9e3cf8e 278
b79bfc3e 279 if(fOldPid){
7ce8802c 280 Double_t mom=track->P();
281 AliAODPid *pidObj = track->GetDetPid();
282
283 Double_t dedx=pidObj->GetITSsignal();
913d1957 284 UChar_t clumap=track->GetITSClusterMap();
285 Int_t nPointsForPid=0;
286 for(Int_t i=2; i<6; i++){
287 if(clumap&(1<<i)) ++nPointsForPid;
288 }
289
290 Bool_t isSA=kTRUE;
291 if(track->GetStatus() & AliESDtrack::kTPCin) isSA = kFALSE;
292
7ce8802c 293 AliITSPIDResponse itsResponse;
f9e3cf8e 294 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
228f3aba 295 Double_t nsigmaMax=fnSigma[4];
f9e3cf8e 296 for(Int_t ipart=0;ipart<5;ipart++){
7ce8802c 297 AliPID::EParticleType type=AliPID::EParticleType(ipart);
913d1957 298 Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type,nPointsForPid,isSA));
228f3aba 299 if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
7ce8802c 300 pid=ipart;
301 nsigmaMax=nsigma;
302 }
303 }
304 }else{ // asks only for one particle specie
305 AliPID::EParticleType type=AliPID::EParticleType(specie);
306 Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
228f3aba 307 if (nsigma>fnSigma[4]) {
f9e3cf8e 308 pid=-1;
7ce8802c 309 }else{
79f3c225 310 pid=specie;
7ce8802c 311 }
312 }
b79bfc3e 313 }else{ // old pid
314
315 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
316 Double_t nsigmaMax=fnSigma[4];
317 for(Int_t ipart=0;ipart<5;ipart++){
318 AliPID::EParticleType type=AliPID::EParticleType(ipart);
319 Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasITS(track,type));
320 if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
321 pid=ipart;
322 nsigmaMax=nsigma;
323 }
324 }
325 }else{ // asks only for one particle specie
326 AliPID::EParticleType type=AliPID::EParticleType(specie);
327 Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasITS(track,type));
328 if (nsigma>fnSigma[4]) {
329 pid=-1;
330 }else{
331 pid=specie;
332 }
333 }
334 } //new pid
335
7ce8802c 336 return pid;
337}
338//----------------------------
228f3aba 339Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
340// n-sigma cut, TOF PID
f9e3cf8e 341
70e398d3 342 if(!CheckTOFPIDStatus(track)) return 0;
f9e3cf8e 343
7ce8802c 344 Double_t time[AliPID::kSPECIESN];
b257e565 345 Double_t sigmaTOFPid[AliPID::kSPECIES];
7ce8802c 346 AliAODPid *pidObj = track->GetDetPid();
347 pidObj->GetIntegratedTimes(time);
228f3aba 348 Double_t sigTOF=pidObj->GetTOFsignal();
b257e565 349 pidObj->GetTOFpidResolution(sigmaTOFPid);
350
f9e3cf8e 351 Int_t pid=-1;
7ce8802c 352
b257e565 353 if(specie<0){
354 Double_t sigmaTOFtrack;
355 if (sigmaTOFPid[4]>0) sigmaTOFtrack=sigmaTOFPid[4];
356 else sigmaTOFtrack=fTOFSigma;
357 Double_t nsigmaMax=sigmaTOFtrack*fnSigma[3];
f9e3cf8e 358 for(Int_t ipart=0;ipart<5;ipart++){
228f3aba 359 Double_t nsigma=TMath::Abs(sigTOF-time[ipart]);
b257e565 360 if (sigmaTOFPid[ipart]>0) sigmaTOFtrack=sigmaTOFPid[ipart];
361 else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
362 if((nsigma<nsigmaMax) && (nsigma<fnSigma[3]*sigmaTOFtrack)) {
7ce8802c 363 pid=ipart;
364 nsigmaMax=nsigma;
365 }
366 }
367 }else{ // asks only for one particle specie
228f3aba 368 Double_t nsigma=TMath::Abs(sigTOF-time[specie]);
b257e565 369 Double_t sigmaTOFtrack;
370 if (sigmaTOFPid[specie]>0) sigmaTOFtrack=sigmaTOFPid[specie];
371 else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
372 if (nsigma>fnSigma[3]*sigmaTOFtrack) {
373 pid=-1;
374 }else{
375 pid=specie;
376 }
7ce8802c 377 }
378 return pid;
379
380}
381//------------------------------
228f3aba 382void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
383// combined PID stored inside the AOD track
7ce8802c 384
385 const Double_t *pid=track->PID();
386 Float_t max=0.;
387 Int_t k=-1;
388 for (Int_t i=0; i<10; i++) {
389 if (pid[i]>max) {k=i; max=pid[i];}
390 }
391
392 if(k==2) type[0]=kTRUE;
393 if(k==3) type[1]=kTRUE;
394 if(k==4) type[2]=kTRUE;
395
396 return;
397}
398//--------------------
cdac45d3 399void AliAODPidHF::BayesianProbability(AliAODTrack *track,Double_t *pid) const{
228f3aba 400// bayesian PID for single detectors or combined
7ce8802c 401
cdac45d3 402 if(fITS && !fTPC && !fTOF) {BayesianProbabilityITS(track,pid);return;}
403 if(fTPC && !fITS && !fTOF) {BayesianProbabilityTPC(track,pid);return;}
404 if(fTOF && !fITS && !fTPC) {BayesianProbabilityTOF(track,pid);return;}
405
406 Double_t probITS[5]={1.,1.,1.,1.,1.};
407 Double_t probTPC[5]={1.,1.,1.,1.,1.};
408 Double_t probTOF[5]={1.,1.,1.,1.,1.};
409 if(fITS) BayesianProbabilityITS(track,probITS);
410 if(fTPC) BayesianProbabilityTPC(track,probTPC);
411 if(fTOF) BayesianProbabilityTOF(track,probTOF);
7ce8802c 412 Double_t probTot[5]={0.,0.,0.,0.,0.};
413 for(Int_t i=0;i<5;i++){
414 probTot[i]=probITS[i]*probTPC[i]*probTOF[i];
415 }
416 for(Int_t i2=0;i2<5;i2++){
417 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]);
418 }
7ce8802c 419
420 return;
421
422}
423//------------------------------------
228f3aba 424void AliAODPidHF::BayesianProbabilityITS(AliAODTrack *track,Double_t *prob) const{
7ce8802c 425
228f3aba 426// bayesian PID for ITS
7ce8802c 427 AliAODpidUtil pid;
428 Double_t itspid[AliPID::kSPECIES];
429 pid.MakeITSPID(track,itspid);
430 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
369ad90e 431 if(fTOF || fTPC || fTRD){
432 prob[ind]=itspid[ind];
433 }else{
434 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]);
435 }
7ce8802c 436 }
437 return;
438
439}
440//------------------------------------
228f3aba 441void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack *track,Double_t *prob) const{
442// bayesian PID for TPC
7ce8802c 443
444 AliAODpidUtil pid;
445 Double_t tpcpid[AliPID::kSPECIES];
446 pid.MakeTPCPID(track,tpcpid);
447 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
369ad90e 448 if(fTOF || fITS || fTRD){
449 prob[ind]=tpcpid[ind];
7ce8802c 450 }else{
369ad90e 451 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]);
7ce8802c 452 }
369ad90e 453}
7ce8802c 454 return;
455
456}
457//------------------------------------
228f3aba 458void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack *track,Double_t *prob) const{
459// bayesian PID for TOF
7ce8802c 460
461 AliAODpidUtil pid;
462 Double_t tofpid[AliPID::kSPECIES];
228f3aba 463 pid.MakeTOFPID(track,tofpid);
7ce8802c 464 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
369ad90e 465 if(fTPC || fITS || fTRD){
466 prob[ind]=tofpid[ind];
467 }else{
7ce8802c 468 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]);
469 }
369ad90e 470}
7ce8802c 471 return;
472
473}
f9e3cf8e 474//---------------------------------
228f3aba 475void AliAODPidHF::BayesianProbabilityTRD(AliAODTrack *track,Double_t *prob) const{
476// bayesian PID for TRD
f9e3cf8e 477
478 AliAODpidUtil pid;
479 Double_t trdpid[AliPID::kSPECIES];
480 pid.MakeTRDPID(track,trdpid);
481 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
369ad90e 482 if(fTPC || fITS || fTOF){
483 prob[ind]=trdpid[ind];
f9e3cf8e 484 }else{
369ad90e 485 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]);
f9e3cf8e 486 }
369ad90e 487}
f9e3cf8e 488 return;
489
490 }
491//--------------------------------
70e398d3 492Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
493 // ITS PID quality cuts
494 if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
495 UChar_t clumap=track->GetITSClusterMap();
496 Int_t nPointsForPid=0;
497 for(Int_t i=2; i<6; i++){
498 if(clumap&(1<<i)) ++nPointsForPid;
499 }
500 if(nPointsForPid<3) return kFALSE;
501 return kTRUE;
502}
503//--------------------------------
504Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
505 // TPC PID quality cuts
506 if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
507 UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
508 if (nTPCClus<70) return kFALSE;
509 return kTRUE;
510}
511//--------------------------------
512Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
513 // TOF PID quality cuts
514 if ((track->GetStatus()&AliESDtrack::kTOFout )==0) return kFALSE;
515 if ((track->GetStatus()&AliESDtrack::kTIME )==0) return kFALSE;
516 if ((track->GetStatus()&AliESDtrack::kTOFpid )==0) return kFALSE;
517 if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0) return kFALSE;
518 return kTRUE;
519}
520//--------------------------------
521Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
522 // TRD PID quality cuts
523 if ((track->GetStatus()&AliESDtrack::kTRDout )==0) return kFALSE;
524 return kTRUE;
525}
526//--------------------------------
228f3aba 527Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
f9e3cf8e 528
228f3aba 529// Quality cuts on the tracks, detector by detector
f9e3cf8e 530
531 if(detectors.Contains("ITS")){
532 if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
533 UChar_t clumap=track->GetITSClusterMap();
534 Int_t nPointsForPid=0;
535 for(Int_t i=2; i<6; i++){
536 if(clumap&(1<<i)) ++nPointsForPid;
537 }
913d1957 538 if(nPointsForPid<3) return kFALSE;
f9e3cf8e 539 }
540
541 if(detectors.Contains("TPC")){
542 if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
543 UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
544 if (nTPCClus<70) return kFALSE;
545 }
546
547 if(detectors.Contains("TOF")){
548 if ((track->GetStatus()&AliESDtrack::kTOFout )==0) return kFALSE;
549 if ((track->GetStatus()&AliESDtrack::kTIME )==0) return kFALSE;
550 if ((track->GetStatus()&AliESDtrack::kTOFpid )==0) return kFALSE;
b257e565 551 if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0) return kFALSE;
f9e3cf8e 552 }
553
554
555 if(detectors.Contains("TRD")){
556 if ((track->GetStatus()&AliESDtrack::kTRDout )==0) return kFALSE;
375df9ce 557 //UChar_t ntracklets = track->GetTRDntrackletsPID();
558 //if(ntracklets<4) return kFALSE;
f9e3cf8e 559 }
560
561 return kTRUE;
562}
228f3aba 563//--------------------------------------------
564Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
565// TPC nsigma cut PID, different sigmas in different p bins
566
70e398d3 567 if(!CheckTPCPIDStatus(track)) return kFALSE;
228f3aba 568 AliAODPid *pidObj = track->GetDetPid();
569 Double_t mom = pidObj->GetTPCmomentum();
2c9eefd4 570 if(mom>fPtThresholdTPC) return 0;
b79bfc3e 571 Double_t nsigma=999.;
572 if(fOldPid){
70e398d3 573 Double_t dedx=pidObj->GetTPCsignal();
574 UShort_t nTPCClus=pidObj->GetTPCsignalN();
575 if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
228f3aba 576
70e398d3 577 AliPID::EParticleType type=AliPID::EParticleType(specie);
578 nsigma = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
b79bfc3e 579
580 }else{ //old pid
70e398d3 581 AliPID::EParticleType type=AliPID::EParticleType(specie);
582 nsigma = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
583 } //new pid
228f3aba 584
585 if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
586 if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
587 if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
588
70e398d3 589 return kFALSE;
228f3aba 590}
591//------------------
70e398d3 592Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
593 // combination of the PID info coming from TPC and TOF
228f3aba 594
70e398d3 595 Bool_t okTPC=CheckTPCPIDStatus(track);
596 Bool_t okTOF=CheckTOFPIDStatus(track);
228f3aba 597
70e398d3 598 if(fMatch==1){
599 //TOF || TPC (a la' Andrea R.)
600 // convention:
601 // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
602 // the method returns the sum of the response of the 2 detectors
228f3aba 603
70e398d3 604 if(fTPC && fTOF) {
605 if(!okTPC && !okTOF) return 0;
606 }
607
608 Int_t tTPCinfo=0;
609 if(fTPC && okTPC){
610 tTPCinfo=-1;
611 if(fAsym) {
612 if(TPCRawAsym(track,specie)) tTPCinfo=1;
613 }else{
614 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
615 }
616 if(fCompat && tTPCinfo<0){
617 Double_t sig0tmp=fnSigma[0];
618 SetSigma(0,fnSigmaCompat[0]);
619 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
620 SetSigma(0,sig0tmp);
621 }
622 }
228f3aba 623
70e398d3 624 Int_t tTOFinfo=0;
625 if(fTOF){
626 if(!okTOF && fTPC) return tTPCinfo;
627 tTOFinfo=-1;
628 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
629 if(fCompat && tTOFinfo>0){
630 Double_t ptrack=track->P();
631 if(ptrack>fPCompatTOF) {
632 Double_t sig0tmp=fnSigma[3];
633 SetSigma(3,fnSigmaCompat[1]);
634 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=0;
635 SetSigma(3,sig0tmp);
636 }
637 }
638 }
913d1957 639
439d61b7 640
70e398d3 641 if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
642 if(!okTOF) return tTPCinfo;
643 return tTOFinfo;
644 }
645
646 if(tTPCinfo+tTOFinfo==0 && fITS){
647 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
648 Int_t tITSinfo = -1;
649 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
650 return tITSinfo;
651 }
652 return tTPCinfo+tTOFinfo;
cdac45d3 653 }
913d1957 654
70e398d3 655 if(fMatch==2){
656 //TPC & TOF (a la' Yifei)
657 // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
658 Int_t tTPCinfo=0;
228f3aba 659
70e398d3 660 if(fTPC && okTPC) {
661 tTPCinfo=1;
662 if(fAsym){
663 if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
664 }else{
665 if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
666 }
667 }
228f3aba 668
70e398d3 669 Int_t tTOFinfo=1;
670 if(fTOF){
671 if(fTPC && !okTOF) return tTPCinfo;
672 if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
673 }
228f3aba 674
70e398d3 675 if(tTOFinfo==1 && tTPCinfo==1) return 1;
228f3aba 676
70e398d3 677 if(tTPCinfo+tTOFinfo==0 && fITS){
678 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
679 Int_t tITSinfo = -1;
680 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
681 return tITSinfo;
682 }
683 return -1;
684 }
913d1957 685
70e398d3 686 if(fMatch==3){
687 //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
688 // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
689 if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
690
691 Double_t ptrack=track->P();
692
693 Int_t tTPCinfo=-1;
694 if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
695 if(!okTPC) return 0;
696 if(fAsym) {
697 if(TPCRawAsym(track,specie)) tTPCinfo=1;
698 }else{
699 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
700 }
701 return tTPCinfo;
702 }
703
704 Int_t tTOFinfo=-1;
705 if(ptrack>=fPLimit[1] && fTOF){
706 if(!okTOF) return 0;
707 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
708 return tTOFinfo;
709 }
710
711 Int_t tITSinfo=-1;
712 if(ptrack<fPLimit[0] && fITS){
713 if(!CheckITSPIDStatus(track)) return 0;
714 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
715 return tITSinfo;
716 }
717 }
228f3aba 718
70e398d3 719 return -1;
cdac45d3 720
721}
722//----------------------------------
723Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
1cb0f658 724// general method to compute PID
70e398d3 725 if(fMatch>0){
726 return MatchTPCTOF(track,specie);
727 }else{
728 if(fTPC && !fTOF && !fITS) {
729 Int_t tTPCres=0;
730 if(!fAsym){
731 tTPCres=ApplyPidTPCRaw(track,specie);
732 if(tTPCres==specie) return 1;
733 else return tTPCres;
734 }else{
735 if(TPCRawAsym(track,specie)) tTPCres=1;
736 else tTPCres=-1;
737 }
738 return tTPCres;
739 }else if(fTOF && !fTPC && !fITS) {
740 Int_t tTOFres=ApplyPidTOFRaw(track,specie);
741 if(tTOFres==specie) return 1;
742 else return tTOFres;
743 }else if(fITS && !fTPC && !fTOF) {
744 Int_t tITSres=ApplyPidITSRaw(track,specie);
745 if(tITSres==specie) return 1;
746 else return tITSres;
ad42e35b 747 }else{
70e398d3 748 AliError("You should enable just one detector if you don't want to match");
749 return 0;
ad42e35b 750 }
70e398d3 751 }
228f3aba 752}
12141139 753//--------------------------------------------
70e398d3 754void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
755 // TPC bethe bloch parameters
32a978c2 756 if(fMC) { // MC
12141139 757
32a978c2 758 if(fPbPb) { // PbPb MC
759
9ea811e0 760 alephParameters[0] = 1.44405/50.;
32a978c2 761 alephParameters[1] = 2.35409e+01;
762 alephParameters[2] = TMath::Exp(-2.90330e+01);
763 alephParameters[3] = 2.10681e+00;
764 alephParameters[4] = 4.62254e+00;
12141139 765
32a978c2 766 } else { // pp MC
5821ea5b 767 if(fMCLowEn2011){
768 alephParameters[0]=0.0207667;
769 alephParameters[1]=29.9936;
770 alephParameters[2]=3.87866e-11;
771 alephParameters[3]=2.17291;
772 alephParameters[4]=7.1623;
3d4afb8c 773 }else if(fOnePad){
774 alephParameters[0]=0.029021;
775 alephParameters[1]=25.4181;
776 alephParameters[2]=4.66596e-08;
777 alephParameters[3]=1.90008;
778 alephParameters[4]=4.63783;
5821ea5b 779 }else{
780 alephParameters[0] = 2.15898/50.;
781 alephParameters[1] = 1.75295e+01;
782 alephParameters[2] = 3.40030e-09;
783 alephParameters[3] = 1.96178e+00;
784 alephParameters[4] = 3.91720e+00;
785 }
32a978c2 786 }
12141139 787
32a978c2 788 } else { // Real Data
789
790 if(fOnePad) { // pp 1-pad (since LHC10d)
791
67a82ef2 792 alephParameters[0] =1.34490e+00/50.;
793 alephParameters[1] = 2.69455e+01;
794 alephParameters[2] = TMath::Exp(-2.97552e+01);
795 alephParameters[3] = 2.35339e+00;
796 alephParameters[4] = 5.98079e+00;
797
798 } else if(fPbPb) { // PbPb
799
800 // alephParameters[0] = 1.25202/50.;
801 // alephParameters[1] = 2.74992e+01;
802 // alephParameters[2] = TMath::Exp(-3.31517e+01);
803 // alephParameters[3] = 2.46246;
804 // alephParameters[4] = 6.78938;
805
806 alephParameters[0] = 5.10207e+00/50.;
807 alephParameters[1] = 7.94982e+00;
808 alephParameters[2] = TMath::Exp(-9.07942e+00);
809 alephParameters[3] = 2.38808e+00;
810 alephParameters[4] = 1.68165e+00;
811
812 } else if(fppLowEn2011){ // pp low energy
813
814 alephParameters[0]=0.031642;
815 alephParameters[1]=22.353;
816 alephParameters[2]=4.16239e-12;
817 alephParameters[3]=2.61952;
818 alephParameters[4]=5.76086;
819
820 } else { // pp no 1-pad (LHC10bc)
821
822 alephParameters[0] = 0.0283086/0.97;
823 alephParameters[1] = 2.63394e+01;
824 alephParameters[2] = 5.04114e-11;
825 alephParameters[3] = 2.12543e+00;
826 alephParameters[4] = 4.88663e+00;
827
828 }
98db0ef2 829
12141139 830 }
831
70e398d3 832}
833
834//-----------------------
835void AliAODPidHF::SetBetheBloch() {
836 // Set Bethe Bloch Parameters
12141139 837
70e398d3 838 Double_t alephParameters[5];
839 GetTPCBetheBlochParams(alephParameters);
840 fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
12141139 841
842 return;
12141139 843}
d50b25df 844//-----------------------
845Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
846
847
70e398d3 848 if(!CheckTOFPIDStatus(track)) return 0;
d50b25df 849
850 Double_t time[AliPID::kSPECIESN];
851 Double_t sigmaTOFPid[AliPID::kSPECIES];
852 AliAODPid *pidObj = track->GetDetPid();
853 pidObj->GetIntegratedTimes(time);
854 Double_t sigTOF=pidObj->GetTOFsignal();
855 pidObj->GetTOFpidResolution(sigmaTOFPid);
856 Double_t sigmaTOFtrack;
857 if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
858 else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
859
860 if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
861
862 return kFALSE;
863
864}
64e39959 865//--------------------------------------------------------------------------
866void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
867
868 //
869 // method setting the prior distributions to the AliPIDCombined object of the AliAODPidHF data member
870 // all the checks are done directly in the AliPIDCombined object
871 //
872
873 GetPidCombined()->SetPriorDistribution(type,prior);
874}
875//--------------------------------------------------------------------------
876void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
877
878 //
879 // Drawing prior distribution for type "type"
880
881 new TCanvas();
882 GetPidCombined()->GetPriorDistribution(type)->Draw();
883}
d50b25df 884
451c9a4e 885//--------------------------------------------------------------------------
886Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const{
887
70e398d3 888 if(!CheckTPCPIDStatus(track)) return -1;
451c9a4e 889
890 Double_t nsigmaTPC=-999;
891
892 if(fOldPid){
893 AliAODPid *pidObj = track->GetDetPid();
894 Double_t dedx=pidObj->GetTPCsignal();
895 Double_t mom = pidObj->GetTPCmomentum();
896 if(mom>fPtThresholdTPC) return -2;
451c9a4e 897 UShort_t nTPCClus=pidObj->GetTPCsignalN();
898 if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
899 AliPID::EParticleType type=AliPID::EParticleType(species);
70e398d3 900 nsigmaTPC = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
451c9a4e 901 sigma=nsigmaTPC;
902 } else{
903
904 AliPID::EParticleType type=AliPID::EParticleType(species);
905 nsigmaTPC = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
906 sigma=nsigmaTPC;
907 }
908 return 1;
909}
910
911//-----------------------------
912
913Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &sigma) const{
914
70e398d3 915 if(!CheckTOFPIDStatus(track)) return -1;
451c9a4e 916
917 Double_t time[AliPID::kSPECIESN];
918 Double_t sigmaTOFPid[AliPID::kSPECIES];
919 AliAODPid *pidObj = track->GetDetPid();
920 pidObj->GetIntegratedTimes(time);
921 Double_t sigTOF=pidObj->GetTOFsignal();
922 pidObj->GetTOFpidResolution(sigmaTOFPid);
923
924 if(sigmaTOFPid[species]<1e-99) return -2;
925
926 Double_t sigmaTOF=TMath::Abs(sigTOF-time[species])/sigmaTOFPid[species];
927 sigma=sigmaTOF;
928 return 1;
929}
930
931//-----------------------------