]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAODPidHF.cxx
Merge branch master into TRDdev
[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>
4c177f98 25#include <TString.h>
26#include <TH1F.h>
27#include <TF1.h>
28#include <TFile.h>
64e39959 29
7ce8802c 30#include "AliAODPidHF.h"
31#include "AliAODPid.h"
7ce8802c 32#include "AliPID.h"
b79bfc3e 33#include "AliPIDResponse.h"
7ce8802c 34#include "AliAODpidUtil.h"
f9e3cf8e 35#include "AliESDtrack.h"
7ce8802c 36
37
38ClassImp(AliAODPidHF)
39
40//------------------------------
41AliAODPidHF::AliAODPidHF():
fcff8ce7 42 TObject(),
228f3aba 43 fnNSigma(5),
0add1d8c 44 fnSigma(0),
228f3aba 45 fTOFSigma(160.),
fcff8ce7 46 fCutTOFmismatch(0.01),
47 fMinNClustersTPCPID(0),
228f3aba 48 fnPriors(5),
0add1d8c 49 fPriors(0),
228f3aba 50 fnPLimit(2),
0add1d8c 51 fPLimit(0),
cdac45d3 52 fAsym(kFALSE),
53 fTPC(kFALSE),
54 fTOF(kFALSE),
55 fITS(kFALSE),
56 fTRD(kFALSE),
79f3c225 57 fMatch(0),
913d1957 58 fCompat(kFALSE),
ad42e35b 59 fPCompatTOF(1.5),
563bd085 60 fUseAsymTOF(kFALSE),
61 fLownSigmaTOF(-3.),
62 fUpnSigmaTOF(3.),
63 fLownSigmaCompatTOF(-3.),
64 fUpnSigmaCompatTOF(3.),
02baac36 65 fnNSigmaCompat(2),
0add1d8c 66 fnSigmaCompat(0),
12141139 67 fMC(kFALSE),
17c4a8d5 68 fOnePad(kFALSE),
5821ea5b 69 fMCLowEn2011(kFALSE),
67a82ef2 70 fppLowEn2011(kFALSE),
439d61b7 71 fPbPb(kFALSE),
b79bfc3e 72 fTOFdecide(kFALSE),
0cda1541 73 fOldPid(kFALSE),
2c9eefd4 74 fPtThresholdTPC(999999.),
c5541cc2 75 fMaxTrackMomForCombinedPID(999999.),
11690a06 76 fPidResponse(0),
70e398d3 77 fPidCombined(new AliPIDCombined()),
4c177f98 78 fTPCResponse(new AliTPCPIDResponse()),
79 fPriorsH(),
80 fCombDetectors(kTPCTOF),
e51f1625 81 fUseCombined(kFALSE),
82 fDefaultPriors(kTRUE)
7ce8802c 83{
fcff8ce7 84 //
85 // Default constructor
86 //
87 fPLimit=new Double_t[fnPLimit];
88 fnSigma=new Double_t[fnNSigma];
89 fPriors=new Double_t[fnPriors];
90 fnSigmaCompat=new Double_t[fnNSigmaCompat];
91
92 for(Int_t i=0;i<fnNSigma;i++){
93 fnSigma[i]=0.;
94 }
95 for(Int_t i=0;i<fnPriors;i++){
96 fPriors[i]=0.;
97 }
98 for(Int_t i=0;i<fnPLimit;i++){
99 fPLimit[i]=0.;
100 }
101 for(Int_t i=0;i<fnNSigmaCompat;i++){
102 fnSigmaCompat[i]=3.;
103 }
104 for(Int_t i=0; i<3; i++){ // pi, K, proton
105 fMaxnSigmaCombined[i]=3.;
106 fMinnSigmaTPC[i]=-3;
107 fMaxnSigmaTPC[i]=3;
108 fMinnSigmaTOF[i]=-3;
109 fMaxnSigmaTOF[i]=3;
110 }
7ce8802c 111
112}
113//----------------------
114AliAODPidHF::~AliAODPidHF()
115{
fa52905d 116 // destructor
a7d74abb 117 if(fPLimit) delete [] fPLimit;
118 if(fnSigma) delete [] fnSigma;
119 if(fPriors) delete [] fPriors;
120 if(fnSigmaCompat) delete [] fnSigmaCompat;
fa52905d 121 delete fPidCombined;
122
70e398d3 123 delete fTPCResponse;
4c177f98 124 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
125 delete fPriorsH[ispecies];
126 }
7ce8802c 127}
128//------------------------
129AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
fcff8ce7 130 TObject(),
228f3aba 131 fnNSigma(pid.fnNSigma),
375df9ce 132 fnSigma(0),
228f3aba 133 fTOFSigma(pid.fTOFSigma),
fcff8ce7 134 fCutTOFmismatch(pid.fCutTOFmismatch),
135 fMinNClustersTPCPID(pid.fMinNClustersTPCPID),
228f3aba 136 fnPriors(pid.fnPriors),
375df9ce 137 fPriors(0),
228f3aba 138 fnPLimit(pid.fnPLimit),
375df9ce 139 fPLimit(0),
cdac45d3 140 fAsym(pid.fAsym),
141 fTPC(pid.fTPC),
142 fTOF(pid.fTOF),
143 fITS(pid.fITS),
144 fTRD(pid.fTRD),
145 fMatch(pid.fMatch),
913d1957 146 fCompat(pid.fCompat),
ad42e35b 147 fPCompatTOF(pid.fPCompatTOF),
563bd085 148 fUseAsymTOF(pid.fUseAsymTOF),
149 fLownSigmaTOF(pid.fLownSigmaTOF),
150 fUpnSigmaTOF(pid.fUpnSigmaTOF),
151 fLownSigmaCompatTOF(pid.fLownSigmaCompatTOF),
152 fUpnSigmaCompatTOF(pid.fUpnSigmaCompatTOF),
02baac36 153 fnNSigmaCompat(pid.fnNSigmaCompat),
2d070486 154 fnSigmaCompat(0x0),
12141139 155 fMC(pid.fMC),
17c4a8d5 156 fOnePad(pid.fOnePad),
5821ea5b 157 fMCLowEn2011(pid.fMCLowEn2011),
67a82ef2 158 fppLowEn2011(pid.fppLowEn2011),
439d61b7 159 fPbPb(pid.fPbPb),
b79bfc3e 160 fTOFdecide(pid.fTOFdecide),
161 fOldPid(pid.fOldPid),
2c9eefd4 162 fPtThresholdTPC(pid.fPtThresholdTPC),
c5541cc2 163 fMaxTrackMomForCombinedPID(pid.fMaxTrackMomForCombinedPID),
2d070486 164 fPidResponse(0x0),
165 fPidCombined(0x0),
4c177f98 166 fTPCResponse(0x0),
167 fCombDetectors(pid.fCombDetectors),
e51f1625 168 fUseCombined(pid.fUseCombined),
169 fDefaultPriors(pid.fDefaultPriors)
375df9ce 170{
7ce8802c 171
2d070486 172 fnSigmaCompat=new Double_t[fnNSigmaCompat];
173 for(Int_t i=0;i<fnNSigmaCompat;i++){
174 fnSigmaCompat[i]=pid.fnSigmaCompat[i];
175 }
375df9ce 176 fnSigma = new Double_t[fnNSigma];
177 for(Int_t i=0;i<fnNSigma;i++){
178 fnSigma[i]=pid.fnSigma[i];
179 }
180 fPriors = new Double_t[fnPriors];
181 for(Int_t i=0;i<fnPriors;i++){
7ce8802c 182 fPriors[i]=pid.fPriors[i];
183 }
375df9ce 184 fPLimit = new Double_t[fnPLimit];
185 for(Int_t i=0;i<fnPLimit;i++){
186 fPLimit[i]=pid.fPLimit[i];
7ce8802c 187 }
4c177f98 188 fPriors = new Double_t[fnPriors];
189 for(Int_t i=0;i<fnPriors;i++){
190 fPriors[i]=pid.fPriors[i];
191 }
192 for(Int_t i=0;i<AliPID::kSPECIES;i++){
193 fPriorsH[i]=pid.fPriorsH[i];
194 }
fcff8ce7 195 for(Int_t i=0; i<3; i++){ // pi, K, proton
196 fMaxnSigmaCombined[i]=pid.fMaxnSigmaCombined[i];
197 fMinnSigmaTPC[i]=pid.fMinnSigmaTPC[i];
198 fMaxnSigmaTPC[i]=pid.fMaxnSigmaTPC[i];
199 fMinnSigmaTOF[i]=pid.fMinnSigmaTOF[i];
200 fMaxnSigmaTOF[i]=pid.fMaxnSigmaTOF[i];
201 }
7ce8802c 202
2d070486 203 // if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
204 fTPCResponse = new AliTPCPIDResponse();
205 SetBetheBloch();
206 fPidCombined = new AliPIDCombined();
375df9ce 207 //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
208 //fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
fcff8ce7 209
375df9ce 210}
7ce8802c 211//----------------------
228f3aba 212Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
fcff8ce7 213 // raw PID for single detectors, returns the particle type with smaller sigma
214 Int_t specie=-1;
215 if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
216 if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
217 if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
218
f9e3cf8e 219 return specie;
7ce8802c 220
221}
222//---------------------------
228f3aba 223Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
fcff8ce7 224 // checks if the track can be a kaon, raw PID applied for single detectors
225 Int_t specie=0;
7ce8802c 226
fcff8ce7 227 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
228 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
229 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
230
231 if(specie==3) return kTRUE;
232 return kFALSE;
7ce8802c 233}
234//---------------------------
228f3aba 235Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
fcff8ce7 236 // checks if the track can be a pion, raw PID applied for single detectors
7ce8802c 237
fcff8ce7 238 Int_t specie=0;
239
240 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
241 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
242 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
243
244 if(specie==2) return kTRUE;
245 return kFALSE;
7ce8802c 246}
247//---------------------------
228f3aba 248Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
fcff8ce7 249 // checks if the track can be a proton raw PID applied for single detectors
250
251 Int_t specie=0;
252 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
253 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
254 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
255
256 if(specie==4) return kTRUE;
257
258 return kFALSE;
7ce8802c 259}
f9e3cf8e 260//--------------------------
228f3aba 261Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
fcff8ce7 262 // checks if the track can be an electron raw PID applied for single detectors
f9e3cf8e 263
fcff8ce7 264 Int_t specie=-1;
265 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
266 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
267 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
268
269 if(specie==0) return kTRUE;
270
271 return kFALSE;
f9e3cf8e 272}
273//--------------------------
228f3aba 274Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
275// n-sigma cut, TPC PID
7ce8802c 276
fcff8ce7 277 Double_t nsigma=-999.;
b79bfc3e 278 Int_t pid=-1;
12141139 279
f9e3cf8e 280 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
4940d5bf 281 Double_t nsigmaMin=999.;
282 for(Int_t ipart=0;ipart<5;ipart++){
283 if(GetnSigmaTPC(track,ipart,nsigma)==1){
284 nsigma=TMath::Abs(nsigma);
285 if((nsigma<nsigmaMin) && (nsigma<fnSigma[0])) {
286 pid=ipart;
287 nsigmaMin=nsigma;
288 }
289 }
7ce8802c 290 }
7ce8802c 291 }else{ // asks only for one particle specie
4940d5bf 292 if(GetnSigmaTPC(track,specie,nsigma)==1){
293 nsigma=TMath::Abs(nsigma);
294 if (nsigma>fnSigma[0]) pid=-1;
295 else pid=specie;
b79bfc3e 296 }
b79bfc3e 297 }
298
4940d5bf 299 return pid;
7ce8802c 300}
301//----------------------------
228f3aba 302Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
303// truncated mean, ITS PID
fcff8ce7 304
305 Double_t nsigma=-999.;
b79bfc3e 306 Int_t pid=-1;
fcff8ce7 307
f9e3cf8e 308 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
fcff8ce7 309 Double_t nsigmaMin=999.;
310 for(Int_t ipart=0;ipart<5;ipart++){
311 if(GetnSigmaITS(track,ipart,nsigma)==1){
312 nsigma=TMath::Abs(nsigma);
313 if((nsigma<nsigmaMin) && (nsigma<fnSigma[4])) {
314 pid=ipart;
315 nsigmaMin=nsigma;
316 }
317 }
7ce8802c 318 }
7ce8802c 319 }else{ // asks only for one particle specie
fcff8ce7 320 if(GetnSigmaITS(track,specie,nsigma)==1){
321 nsigma=TMath::Abs(nsigma);
322 if (nsigma>fnSigma[4]) pid=-1;
323 else pid=specie;
b79bfc3e 324 }
b79bfc3e 325 }
fcff8ce7 326
327 return pid;
7ce8802c 328}
329//----------------------------
228f3aba 330Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
331// n-sigma cut, TOF PID
f9e3cf8e 332
fcff8ce7 333 Double_t nsigma=-999.;
4940d5bf 334 Int_t pid=-1;
f9e3cf8e 335
4940d5bf 336 if(specie<0){
337 Double_t nsigmaMin=999.;
338 for(Int_t ipart=0;ipart<5;ipart++){
339 if(GetnSigmaTOF(track,ipart,nsigma)==1){
340 nsigma=TMath::Abs(nsigma);
341 if((nsigma<nsigmaMin)&& (nsigma<fnSigma[3])){
342 pid=ipart;
343 nsigmaMin=nsigma;
344 }
345 }
346 }
347 }else{ // asks only for one particle specie
563bd085 348 Double_t nSigmaMin,nSigmaMax;
349 if(fUseAsymTOF){
350 nSigmaMin=fLownSigmaTOF;
351 nSigmaMax=fUpnSigmaTOF;
352 }else{
353 nSigmaMin=-fnSigma[3];
354 nSigmaMax=fnSigma[3];
355 }
4940d5bf 356 if(GetnSigmaTOF(track,specie,nsigma)==1){
563bd085 357 if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
4940d5bf 358 else pid=specie;
359 }
360 }
361 return pid;
563bd085 362}
363//----------------------------
364Int_t AliAODPidHF::ApplyTOFCompatibilityBand(AliAODTrack *track,Int_t specie) const{
365// n-sigma cut, TOF PID
7ce8802c 366
563bd085 367 if(specie<0) return -1;
368 Double_t nsigma=-999.;
369 Int_t pid=-1;
370
371 Double_t nSigmaMin,nSigmaMax;
372 if(fUseAsymTOF){
373 nSigmaMin=fLownSigmaCompatTOF;
374 nSigmaMax=fUpnSigmaCompatTOF;
375 }else{
376 nSigmaMin=-fnSigmaCompat[1];
377 nSigmaMax=fnSigmaCompat[1];
7ce8802c 378 }
563bd085 379 if(GetnSigmaTOF(track,specie,nsigma)==1){
380 if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
381 else pid=specie;
382 }
383 return pid;
7ce8802c 384}
385//------------------------------
228f3aba 386void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
387// combined PID stored inside the AOD track
7ce8802c 388
389 const Double_t *pid=track->PID();
390 Float_t max=0.;
391 Int_t k=-1;
392 for (Int_t i=0; i<10; i++) {
393 if (pid[i]>max) {k=i; max=pid[i];}
394 }
395
396 if(k==2) type[0]=kTRUE;
397 if(k==3) type[1]=kTRUE;
398 if(k==4) type[2]=kTRUE;
399
400 return;
401}
f9e3cf8e 402//--------------------------------
70e398d3 403Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
fcff8ce7 404 // Check if the track is good for ITS PID
405 AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kITS,track);
406 if (status != AliPIDResponse::kDetPidOk) return kFALSE;
70e398d3 407 return kTRUE;
408}
409//--------------------------------
410Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
fcff8ce7 411 // Check if the track is good for TPC PID
412 AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTPC,track);
413 if (status != AliPIDResponse::kDetPidOk) return kFALSE;
414 UInt_t nclsTPCPID = track->GetTPCsignalN();
415 if(nclsTPCPID<fMinNClustersTPCPID) return kFALSE;
70e398d3 416 return kTRUE;
417}
418//--------------------------------
419Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
fcff8ce7 420 // Check if the track is good for TOF PID
421 AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTOF,track);
422 if (status != AliPIDResponse::kDetPidOk) return kFALSE;
423 Float_t probMis = fPidResponse->GetTOFMismatchProbability(track);
424 if (probMis > fCutTOFmismatch) return kFALSE;
425 if ((track->GetStatus()&AliESDtrack::kTOFpid )==0 &&
426 track->GetStatus()&AliESDtrack::kITSrefit ) return kFALSE;
70e398d3 427 return kTRUE;
428}
429//--------------------------------
430Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
fcff8ce7 431 // Check if the track is good for TRD PID
432 AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTRD,track);
433 if (status != AliPIDResponse::kDetPidOk) return kFALSE;
70e398d3 434 return kTRUE;
435}
436//--------------------------------
228f3aba 437Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
fcff8ce7 438 // Quality cuts on the tracks, detector by detector
439 if(detectors.Contains("ITS")) return CheckITSPIDStatus(track);
440 else if(detectors.Contains("TPC")) return CheckTPCPIDStatus(track);
441 else if(detectors.Contains("TOF")) return CheckTOFPIDStatus(track);
442 else if(detectors.Contains("TRD")) return CheckTRDPIDStatus(track);
443 else{
444 AliError("Wrong detector name");
445 return kFALSE;
f9e3cf8e 446 }
f9e3cf8e 447}
228f3aba 448//--------------------------------------------
449Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
fcff8ce7 450 // TPC nsigma cut PID, different sigmas in different p bins
228f3aba 451
c5541cc2 452 AliAODPid *pidObj = track->GetDetPid();
453 Double_t mom = pidObj->GetTPCmomentum();
454 if(mom>fPtThresholdTPC) return kTRUE;
455
4940d5bf 456 Double_t nsigma;
457 if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
458 nsigma=TMath::Abs(nsigma);
459
228f3aba 460
461 if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
462 if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
463 if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
464
70e398d3 465 return kFALSE;
228f3aba 466}
467//------------------
70e398d3 468Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
469 // combination of the PID info coming from TPC and TOF
228f3aba 470
c5541cc2 471 Double_t ptrack=track->P();
472 if(ptrack>fMaxTrackMomForCombinedPID) return 1;
473
70e398d3 474 Bool_t okTPC=CheckTPCPIDStatus(track);
c5541cc2 475 if(ptrack>fPtThresholdTPC) okTPC=kFALSE;
fcff8ce7 476 Bool_t okTOF=CheckTOFPIDStatus(track);
228f3aba 477
70e398d3 478 if(fMatch==1){
479 //TOF || TPC (a la' Andrea R.)
480 // convention:
481 // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
482 // the method returns the sum of the response of the 2 detectors
228f3aba 483
70e398d3 484 if(fTPC && fTOF) {
485 if(!okTPC && !okTOF) return 0;
486 }
487
488 Int_t tTPCinfo=0;
489 if(fTPC && okTPC){
490 tTPCinfo=-1;
491 if(fAsym) {
492 if(TPCRawAsym(track,specie)) tTPCinfo=1;
493 }else{
494 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
495 }
496 if(fCompat && tTPCinfo<0){
497 Double_t sig0tmp=fnSigma[0];
498 SetSigma(0,fnSigmaCompat[0]);
499 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
500 SetSigma(0,sig0tmp);
501 }
502 }
228f3aba 503
70e398d3 504 Int_t tTOFinfo=0;
505 if(fTOF){
506 if(!okTOF && fTPC) return tTPCinfo;
507 tTOFinfo=-1;
508 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
509 if(fCompat && tTOFinfo>0){
70e398d3 510 if(ptrack>fPCompatTOF) {
563bd085 511 if(ApplyTOFCompatibilityBand(track,specie)==specie) tTOFinfo=0;
70e398d3 512 }
513 }
514 }
913d1957 515
439d61b7 516
70e398d3 517 if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
518 if(!okTOF) return tTPCinfo;
519 return tTOFinfo;
520 }
521
522 if(tTPCinfo+tTOFinfo==0 && fITS){
523 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
524 Int_t tITSinfo = -1;
525 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
526 return tITSinfo;
527 }
528 return tTPCinfo+tTOFinfo;
cdac45d3 529 }
913d1957 530
70e398d3 531 if(fMatch==2){
532 //TPC & TOF (a la' Yifei)
533 // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
534 Int_t tTPCinfo=0;
228f3aba 535
70e398d3 536 if(fTPC && okTPC) {
537 tTPCinfo=1;
538 if(fAsym){
539 if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
540 }else{
541 if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
542 }
543 }
228f3aba 544
70e398d3 545 Int_t tTOFinfo=1;
546 if(fTOF){
547 if(fTPC && !okTOF) return tTPCinfo;
548 if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
549 }
228f3aba 550
70e398d3 551 if(tTOFinfo==1 && tTPCinfo==1) return 1;
228f3aba 552
70e398d3 553 if(tTPCinfo+tTOFinfo==0 && fITS){
554 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
555 Int_t tITSinfo = -1;
556 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
557 return tITSinfo;
558 }
559 return -1;
560 }
913d1957 561
70e398d3 562 if(fMatch==3){
563 //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
564 // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
565 if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
566
c5541cc2 567
70e398d3 568 Int_t tTPCinfo=-1;
569 if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
570 if(!okTPC) return 0;
571 if(fAsym) {
572 if(TPCRawAsym(track,specie)) tTPCinfo=1;
573 }else{
574 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
575 }
576 return tTPCinfo;
577 }
578
579 Int_t tTOFinfo=-1;
580 if(ptrack>=fPLimit[1] && fTOF){
581 if(!okTOF) return 0;
582 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
583 return tTOFinfo;
584 }
585
586 Int_t tITSinfo=-1;
587 if(ptrack<fPLimit[0] && fITS){
588 if(!CheckITSPIDStatus(track)) return 0;
589 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
590 return tITSinfo;
591 }
592 }
228f3aba 593
fcff8ce7 594 if(fMatch==4 || fMatch==5){
595
596 // fMatch == 4 ---> "circular cut" in nSigmaTPC, nSimgaTOF plane
597 // ---> nsigmaTPC^2+nsigmaTOF^2 < cut^2
598 // fMatch == 5 ---> "rectangular cut" in nSigmaTPC, nsigmaTOF plane
599 // ---> ns1<nSigmaTPC<NS1 && ns2<nSigmaTOF<NS2
600
601 Double_t nSigmaTPC=0.;
602 if(okTPC) {
603 nSigmaTPC=fPidResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)specie);
604 if(nSigmaTPC<-990.) nSigmaTPC=0.;
605 }
606 Double_t nSigmaTOF=0.;
607 if(okTOF) {
608 nSigmaTOF=fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)specie);
609 }
610 Int_t iPart=specie-2; //species is 2 for pions,3 for kaons and 4 for protons
611 if(iPart<0 || iPart>2) return -1;
612 if(fMatch==4){
613 Double_t nSigma2=nSigmaTPC*nSigmaTPC+nSigmaTOF*nSigmaTOF;
614 if(nSigma2<fMaxnSigmaCombined[iPart]*fMaxnSigmaCombined[iPart]) return 1;
615 else return -1;
616 }
617 else if(fMatch==5){
618 if((nSigmaTPC>fMinnSigmaTPC[iPart] && nSigmaTPC<fMaxnSigmaTPC[iPart]) &&
619 (nSigmaTOF>fMinnSigmaTOF[iPart] && nSigmaTOF<fMaxnSigmaTOF[iPart])) return 1;
620 else return -1;
621 }
622 }
623
624
70e398d3 625 return -1;
cdac45d3 626
627}
628//----------------------------------
629Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
fcff8ce7 630 // general method to compute PID
70e398d3 631 if(fMatch>0){
632 return MatchTPCTOF(track,specie);
633 }else{
634 if(fTPC && !fTOF && !fITS) {
635 Int_t tTPCres=0;
636 if(!fAsym){
637 tTPCres=ApplyPidTPCRaw(track,specie);
638 if(tTPCres==specie) return 1;
639 else return tTPCres;
640 }else{
641 if(TPCRawAsym(track,specie)) tTPCres=1;
642 else tTPCres=-1;
643 }
644 return tTPCres;
645 }else if(fTOF && !fTPC && !fITS) {
646 Int_t tTOFres=ApplyPidTOFRaw(track,specie);
647 if(tTOFres==specie) return 1;
648 else return tTOFres;
649 }else if(fITS && !fTPC && !fTOF) {
650 Int_t tITSres=ApplyPidITSRaw(track,specie);
651 if(tITSres==specie) return 1;
652 else return tITSres;
ad42e35b 653 }else{
70e398d3 654 AliError("You should enable just one detector if you don't want to match");
655 return 0;
ad42e35b 656 }
70e398d3 657 }
228f3aba 658}
12141139 659//--------------------------------------------
70e398d3 660void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
661 // TPC bethe bloch parameters
fcff8ce7 662 if(fMC) { // MC
663
664 if(fPbPb) { // PbPb MC
665
666 alephParameters[0] = 1.44405/50.;
667 alephParameters[1] = 2.35409e+01;
668 alephParameters[2] = TMath::Exp(-2.90330e+01);
669 alephParameters[3] = 2.10681e+00;
670 alephParameters[4] = 4.62254e+00;
671
672 } else { // pp MC
673 if(fMCLowEn2011){
674 alephParameters[0]=0.0207667;
675 alephParameters[1]=29.9936;
676 alephParameters[2]=3.87866e-11;
677 alephParameters[3]=2.17291;
678 alephParameters[4]=7.1623;
679 }else if(fOnePad){
680 alephParameters[0]=0.029021;
681 alephParameters[1]=25.4181;
682 alephParameters[2]=4.66596e-08;
683 alephParameters[3]=1.90008;
684 alephParameters[4]=4.63783;
685 }else{
686 alephParameters[0] = 2.15898/50.;
687 alephParameters[1] = 1.75295e+01;
688 alephParameters[2] = 3.40030e-09;
689 alephParameters[3] = 1.96178e+00;
690 alephParameters[4] = 3.91720e+00;
691 }
692 }
693
694 } else { // Real Data
695
696 if(fOnePad) { // pp 1-pad (since LHC10d)
697
698 alephParameters[0] =1.34490e+00/50.;
699 alephParameters[1] = 2.69455e+01;
700 alephParameters[2] = TMath::Exp(-2.97552e+01);
701 alephParameters[3] = 2.35339e+00;
702 alephParameters[4] = 5.98079e+00;
703
704 } else if(fPbPb) { // PbPb
705
706 // alephParameters[0] = 1.25202/50.;
707 // alephParameters[1] = 2.74992e+01;
708 // alephParameters[2] = TMath::Exp(-3.31517e+01);
709 // alephParameters[3] = 2.46246;
710 // alephParameters[4] = 6.78938;
711
712 alephParameters[0] = 5.10207e+00/50.;
713 alephParameters[1] = 7.94982e+00;
714 alephParameters[2] = TMath::Exp(-9.07942e+00);
715 alephParameters[3] = 2.38808e+00;
716 alephParameters[4] = 1.68165e+00;
717
718 } else if(fppLowEn2011){ // pp low energy
719
720 alephParameters[0]=0.031642;
721 alephParameters[1]=22.353;
722 alephParameters[2]=4.16239e-12;
723 alephParameters[3]=2.61952;
724 alephParameters[4]=5.76086;
725
726 } else { // pp no 1-pad (LHC10bc)
727
728 alephParameters[0] = 0.0283086/0.97;
729 alephParameters[1] = 2.63394e+01;
730 alephParameters[2] = 5.04114e-11;
731 alephParameters[3] = 2.12543e+00;
732 alephParameters[4] = 4.88663e+00;
733
734 }
735
736 }
12141139 737
70e398d3 738}
739
740//-----------------------
741void AliAODPidHF::SetBetheBloch() {
742 // Set Bethe Bloch Parameters
12141139 743
70e398d3 744 Double_t alephParameters[5];
745 GetTPCBetheBlochParams(alephParameters);
746 fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
12141139 747
748 return;
12141139 749}
d50b25df 750
4940d5bf 751
64e39959 752//--------------------------------------------------------------------------
fcff8ce7 753Int_t AliAODPidHF::GetnSigmaITS(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
754 // get n sigma for ITS
64e39959 755
64e39959 756
fcff8ce7 757 if (!CheckITSPIDStatus(track)) return -1;
64e39959 758
fcff8ce7 759 Double_t nsigmaITS=-999;
64e39959 760
fcff8ce7 761 if (fOldPid) {
762 Double_t mom=track->P();
763 AliAODPid *pidObj = track->GetDetPid();
764 Double_t dedx=pidObj->GetITSsignal();
765
766 AliITSPIDResponse itsResponse;
767 AliPID::EParticleType type=AliPID::EParticleType(species);
768 nsigmaITS = itsResponse.GetNumberOfSigmas(mom,dedx,type);
769
770 } // old pid
771 else { // new pid
d50b25df 772
fcff8ce7 773 AliPID::EParticleType type=AliPID::EParticleType(species);
774 nsigmaITS = fPidResponse->NumberOfSigmasITS(track,type);
775
776 } //new pid
777
778 nsigma = nsigmaITS;
779
780 return 1;
781
782}
451c9a4e 783//--------------------------------------------------------------------------
4940d5bf 784Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
4c177f98 785 // get n sigma for TPC
786
70e398d3 787 if(!CheckTPCPIDStatus(track)) return -1;
451c9a4e 788
789 Double_t nsigmaTPC=-999;
790
791 if(fOldPid){
792 AliAODPid *pidObj = track->GetDetPid();
793 Double_t dedx=pidObj->GetTPCsignal();
794 Double_t mom = pidObj->GetTPCmomentum();
795 if(mom>fPtThresholdTPC) return -2;
451c9a4e 796 UShort_t nTPCClus=pidObj->GetTPCsignalN();
797 if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
798 AliPID::EParticleType type=AliPID::EParticleType(species);
4940d5bf 799 nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
800 nsigma=nsigmaTPC;
451c9a4e 801 } else{
4940d5bf 802 if(!fPidResponse) return -1;
451c9a4e 803 AliPID::EParticleType type=AliPID::EParticleType(species);
4940d5bf 804 nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
805 nsigma=nsigmaTPC;
451c9a4e 806 }
807 return 1;
808}
809
810//-----------------------------
811
4940d5bf 812Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
4c177f98 813 // get n sigma for TOF
451c9a4e 814
70e398d3 815 if(!CheckTOFPIDStatus(track)) return -1;
4940d5bf 816
817 if(fPidResponse){
818 nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
a585a726 819 return 1;
4940d5bf 820 }else{
a585a726 821 AliFatal("To use TOF PID you need to attach AliPIDResponseTask");
822 nsigma=-999.;
823 return -1;
4940d5bf 824 }
451c9a4e 825}
826
e7af8919 827//-----------------------
4c177f98 828Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
829 // Exclude a given hypothesis (labelTracks) in detector
e7af8919 830
831 if (detectors.Contains("ITS")) {
832
833 AliInfo("Nothing to be done");
834 /*
835 Double_t nsigma=0.;
836 if (GetnSigmaITS(track,labelTrack,nsigma)==1){
837 if(nsigma>nsigmaCut) return kTRUE;
838 }
839 */
840 return kFALSE;
841
842 } else if (detectors.Contains("TPC")) {
843
844 Double_t nsigma=0.;
845 if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
846 if(nsigma>nsigmaCut) return kTRUE;
847 }
848 return kFALSE;
849
850 } else if (detectors.Contains("TOF")) {
851
e7af8919 852 Double_t nsigma=0.;
853 if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
854 if(nsigma>nsigmaCut) return kTRUE;
855 }
856 return kFALSE;
857
858 }
859 return kFALSE;
860
861}
fcff8ce7 862//-----------------------
863Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
864 // TOF proton compatibility
865
866 if(!CheckTOFPIDStatus(track)) return 0;
867
868 Double_t nsigma;
869 if(GetnSigmaTOF(track,3,nsigma)==1){
870 if(nsigma>nsigmaK) return kTRUE;
871 }
872 return kFALSE;
873 /* Double_t time[AliPID::kSPECIESN];
874 Double_t sigmaTOFPid[AliPID::kSPECIES];
875 AliAODPid *pidObj = track->GetDetPid();
876 pidObj->GetIntegratedTimes(time);
877 Double_t sigTOF=pidObj->GetTOFsignal();
878
879 AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
880 if (event) {
881 AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
882 if (tofH && fPidResponse) {
883 AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
884 sigTOF -= TOFres.GetStartTime(track->P());
885 sigmaTOFPid[3]=TOFres.GetExpectedSigma(track->P(),time[3],AliPID::ParticleMass(3));
886 }
887 else pidObj->GetTOFpidResolution(sigmaTOFPid);
888 } else pidObj->GetTOFpidResolution(sigmaTOFPid);
889 Double_t sigmaTOFtrack;
890 if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
891 else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
892
893 if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
894
895 return kFALSE;
896 */
897}
898
899//--------------------------------------------------------------------------
900void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
901
902 //
903 // method setting the prior distributions to the AliPIDCombined object of the AliAODPidHF data member
904 // all the checks are done directly in the AliPIDCombined object
905 //
906
907 GetPidCombined()->SetPriorDistribution(type,prior);
908}
909//--------------------------------------------------------------------------
910void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
911
912 //
913 // Drawing prior distribution for type "type"
914
915 new TCanvas();
916 GetPidCombined()->GetPriorDistribution(type)->Draw();
917}
918
451c9a4e 919//-----------------------------
4c177f98 920void AliAODPidHF::SetPriorsHistos(TString priorFileName){
921 // Set histograms with priors
03a82f54 922
4c177f98 923 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
924 if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
925 TString nt ="name";
926 nt+="_prior_";
927 nt+=AliPID::ParticleName(ispecies);
4c177f98 928 }
03a82f54 929 TDirectory *current = gDirectory;
4c177f98 930 TFile *priorFile=TFile::Open(priorFileName);
931 if (priorFile) {
03a82f54 932 TH1F* h3=static_cast<TH1F*>(priorFile->Get("priors3step9"));
933 TH1F* h2=static_cast<TH1F*>(priorFile->Get("priors2step9"));
934 TH1F* h1=static_cast<TH1F*>(priorFile->Get("priors1step9"));
935 current->cd();
936 fPriorsH[AliPID::kProton] = new TH1F(*h3);
937 fPriorsH[AliPID::kKaon ] = new TH1F(*h2);
938 fPriorsH[AliPID::kPion ] = new TH1F(*h1);
939 priorFile->Close();
4c177f98 940 delete priorFile;
941 TF1 *salt=new TF1("salt","1.e-10",0,10);
942 fPriorsH[AliPID::kProton]->Add(salt);
943 fPriorsH[AliPID::kKaon ]->Add(salt);
944 fPriorsH[AliPID::kPion ]->Add(salt);
945 delete salt;
946 }
947}
948//----------------------------------
949void AliAODPidHF::SetUpCombinedPID(){
950 // Configuration of combined Bayesian PID
951
fcff8ce7 952 fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
e51f1625 953 if(!fDefaultPriors){
954 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
955 fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
956 }
957 }else{
958 fPidCombined->SetDefaultTPCPriors();
4c177f98 959 }
fcff8ce7 960 switch (fCombDetectors){
4c177f98 961 case kTPCTOF:
fcff8ce7 962 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
963 break;
71391326 964 case kTPCAndTOF:
965 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC & AliPIDResponse::kDetTOF);
966 break;
4c177f98 967 case kTPCITS:
fcff8ce7 968 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
969 break;
4c177f98 970 case kTPC:
fcff8ce7 971 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
972 break;
4c177f98 973 case kTOF:
fcff8ce7 974 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
975 break;
976 }
4c177f98 977}
f35c481f 978
f35c481f 979
f35149be 980//-----------------------------
981void AliAODPidHF::PrintAll() const {
982 // print the configuration
983 printf("Detectors used for PID: ");
984 if(fITS) printf("ITS ");
985 if(fTPC) printf("TPC ");
986 if(fTRD) printf("TRD ");
987 if(fTOF) printf("TOF ");
988 printf("\n");
989 printf("Minimum TPC PID clusters = %d\n",fMinNClustersTPCPID);
990 printf("Maximum momentum for using TPC PID = %f\n",fPtThresholdTPC);
991 printf("TOF Mismatch probablility cut = %f\n",fCutTOFmismatch);
992 printf("Maximum momentum for combined PID TPC PID = %f\n",fMaxTrackMomForCombinedPID);
993 if(fOldPid){
994 printf("Use OLD PID");
995 printf(" fMC = %d\n",fMC);
996 printf(" fPbPb = %d\n",fPbPb);
997 printf(" fOnePad = %d\n",fOnePad);
998 printf(" fMCLowEn2011 = %d\n",fMCLowEn2011);
999 printf(" fppLowEn2011 = %d\n",fppLowEn2011);
1000 }
1001 printf("--- Matching algorithm = %d ---\n",fMatch);
1002 if(fMatch==1){
1003 if(fITS) printf("nSigmaITS = %.2f\n",fnSigma[4]);
1004 if(fTOF){
1005 printf("nSigmaTOF = %.2f\n",fnSigma[3]);
1006 if(fCompat) printf("Compatibility band at nSigmaTOF=%.2f for p>%.2f\n",fnSigmaCompat[1],fPCompatTOF);
1007 }
1008 if(fTPC){
1009 if(fAsym){
1010 printf("nSigmaTPC:\n");
1011 printf(" pt<%.2f \t nsigmaTPC= %.2f\n",fPLimit[0],fnSigma[0]);
1012 printf(" %.2f<pt<%.2f \t nsigmaTPC= %.2f\n",fPLimit[0],fPLimit[1],fnSigma[1]);
1013 printf(" pt>%.2f \t nsigmaTPC= %.2f\n",fPLimit[1],fnSigma[2]);
1014 }else{
1015 printf("nSigmaTPC = %.2f\n",fnSigma[0]);
1016 }
1017 if(fCompat) printf("Compatibility band at nSigmaTPC=%.2f\n",fnSigmaCompat[0]);
1018 }
1019 }else if(fMatch==4){
1020 printf("Cuts on sqrt(nSigmaTPC^2+nSigmaTOF^2):\n");
1021 printf(" Pions: nSigma = %.2f\n",fMaxnSigmaCombined[0]);
1022 printf(" Kaons: nSigma = %.2f\n",fMaxnSigmaCombined[1]);
1023 printf(" Protons: nSigma = %.2f\n",fMaxnSigmaCombined[2]);
1024 }else if(fMatch==5){
1025 printf("nSigma ranges:\n");
1026 printf(" Pions: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1027 fMinnSigmaTPC[0],fMaxnSigmaTPC[0],fMinnSigmaTOF[0],fMaxnSigmaTOF[0]);
1028 printf(" Kaons: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1029 fMinnSigmaTPC[1],fMaxnSigmaTPC[1],fMinnSigmaTOF[1],fMaxnSigmaTOF[1]);
1030 printf(" Protons: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1031 fMinnSigmaTPC[2],fMaxnSigmaTPC[2],fMinnSigmaTOF[2],fMaxnSigmaTOF[2]);
1032 }
1033}