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