1 /**************************************************************************
2 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 * *************************************************************************/
19 //***********************************************************
21 // class for PID with AliAODRecoDecayHF
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
23 //***********************************************************
30 #include "AliAODPidHF.h"
31 #include "AliAODPid.h"
33 #include "AliPIDResponse.h"
34 #include "AliAODpidUtil.h"
35 #include "AliESDtrack.h"
40 //------------------------------
41 AliAODPidHF::AliAODPidHF():
46 fCutTOFmismatch(0.01),
47 fMinNClustersTPCPID(0),
63 fLownSigmaCompatTOF(-3.),
64 fUpnSigmaCompatTOF(3.),
74 fPtThresholdTPC(999999.),
75 fMaxTrackMomForCombinedPID(999999.),
77 fPidCombined(new AliPIDCombined()),
78 fTPCResponse(new AliTPCPIDResponse()),
80 fCombDetectors(kTPCTOF),
85 // Default constructor
87 fPLimit=new Double_t[fnPLimit];
88 fnSigma=new Double_t[fnNSigma];
89 fPriors=new Double_t[fnPriors];
90 fnSigmaCompat=new Double_t[fnNSigmaCompat];
92 for(Int_t i=0;i<fnNSigma;i++){
95 for(Int_t i=0;i<fnPriors;i++){
98 for(Int_t i=0;i<fnPLimit;i++){
101 for(Int_t i=0;i<fnNSigmaCompat;i++){
104 for(Int_t i=0; i<3; i++){ // pi, K, proton
105 fMaxnSigmaCombined[i]=3.;
113 //----------------------
114 AliAODPidHF::~AliAODPidHF()
117 if(fPLimit) delete [] fPLimit;
118 if(fnSigma) delete [] fnSigma;
119 if(fPriors) delete [] fPriors;
120 if(fnSigmaCompat) delete [] fnSigmaCompat;
124 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
125 delete fPriorsH[ispecies];
128 //------------------------
129 AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
131 fnNSigma(pid.fnNSigma),
133 fTOFSigma(pid.fTOFSigma),
134 fCutTOFmismatch(pid.fCutTOFmismatch),
135 fMinNClustersTPCPID(pid.fMinNClustersTPCPID),
136 fnPriors(pid.fnPriors),
138 fnPLimit(pid.fnPLimit),
146 fCompat(pid.fCompat),
147 fPCompatTOF(pid.fPCompatTOF),
148 fUseAsymTOF(pid.fUseAsymTOF),
149 fLownSigmaTOF(pid.fLownSigmaTOF),
150 fUpnSigmaTOF(pid.fUpnSigmaTOF),
151 fLownSigmaCompatTOF(pid.fLownSigmaCompatTOF),
152 fUpnSigmaCompatTOF(pid.fUpnSigmaCompatTOF),
153 fnNSigmaCompat(pid.fnNSigmaCompat),
156 fOnePad(pid.fOnePad),
157 fMCLowEn2011(pid.fMCLowEn2011),
158 fppLowEn2011(pid.fppLowEn2011),
160 fTOFdecide(pid.fTOFdecide),
161 fOldPid(pid.fOldPid),
162 fPtThresholdTPC(pid.fPtThresholdTPC),
163 fMaxTrackMomForCombinedPID(pid.fMaxTrackMomForCombinedPID),
167 fCombDetectors(pid.fCombDetectors),
168 fUseCombined(pid.fUseCombined),
169 fDefaultPriors(pid.fDefaultPriors)
172 fnSigmaCompat=new Double_t[fnNSigmaCompat];
173 for(Int_t i=0;i<fnNSigmaCompat;i++){
174 fnSigmaCompat[i]=pid.fnSigmaCompat[i];
176 fnSigma = new Double_t[fnNSigma];
177 for(Int_t i=0;i<fnNSigma;i++){
178 fnSigma[i]=pid.fnSigma[i];
180 fPriors = new Double_t[fnPriors];
181 for(Int_t i=0;i<fnPriors;i++){
182 fPriors[i]=pid.fPriors[i];
184 fPLimit = new Double_t[fnPLimit];
185 for(Int_t i=0;i<fnPLimit;i++){
186 fPLimit[i]=pid.fPLimit[i];
188 fPriors = new Double_t[fnPriors];
189 for(Int_t i=0;i<fnPriors;i++){
190 fPriors[i]=pid.fPriors[i];
192 for(Int_t i=0;i<AliPID::kSPECIES;i++){
193 fPriorsH[i]=pid.fPriorsH[i];
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];
203 // if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
204 fTPCResponse = new AliTPCPIDResponse();
206 fPidCombined = new AliPIDCombined();
207 //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
208 //fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
211 //----------------------
212 Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
213 // raw PID for single detectors, returns the particle type with smaller sigma
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);
222 //---------------------------
223 Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
224 // checks if the track can be a kaon, raw PID applied for single detectors
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);
231 if(specie==3) return kTRUE;
234 //---------------------------
235 Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
236 // checks if the track can be a pion, raw PID applied for single detectors
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);
244 if(specie==2) return kTRUE;
247 //---------------------------
248 Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
249 // checks if the track can be a proton raw PID applied for single detectors
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);
256 if(specie==4) return kTRUE;
260 //--------------------------
261 Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
262 // checks if the track can be an electron raw PID applied for single detectors
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);
269 if(specie==0) return kTRUE;
273 //--------------------------
274 Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
275 // n-sigma cut, TPC PID
277 Double_t nsigma=-999.;
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
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])) {
291 }else{ // asks only for one particle specie
292 if(GetnSigmaTPC(track,specie,nsigma)==1){
293 nsigma=TMath::Abs(nsigma);
294 if (nsigma>fnSigma[0]) pid=-1;
301 //----------------------------
302 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
303 // truncated mean, ITS PID
305 Double_t nsigma=-999.;
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
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])) {
319 }else{ // asks only for one particle specie
320 if(GetnSigmaITS(track,specie,nsigma)==1){
321 nsigma=TMath::Abs(nsigma);
322 if (nsigma>fnSigma[4]) pid=-1;
329 //----------------------------
330 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
331 // n-sigma cut, TOF PID
333 Double_t nsigma=-999.;
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])){
347 }else{ // asks only for one particle specie
348 Double_t nSigmaMin,nSigmaMax;
350 nSigmaMin=fLownSigmaTOF;
351 nSigmaMax=fUpnSigmaTOF;
353 nSigmaMin=-fnSigma[3];
354 nSigmaMax=fnSigma[3];
356 if(GetnSigmaTOF(track,specie,nsigma)==1){
357 if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
363 //----------------------------
364 Int_t AliAODPidHF::ApplyTOFCompatibilityBand(AliAODTrack *track,Int_t specie) const{
365 // n-sigma cut, TOF PID
367 if(specie<0) return -1;
368 Double_t nsigma=-999.;
371 Double_t nSigmaMin,nSigmaMax;
373 nSigmaMin=fLownSigmaCompatTOF;
374 nSigmaMax=fUpnSigmaCompatTOF;
376 nSigmaMin=-fnSigmaCompat[1];
377 nSigmaMax=fnSigmaCompat[1];
379 if(GetnSigmaTOF(track,specie,nsigma)==1){
380 if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
385 //------------------------------
386 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
387 // combined PID stored inside the AOD track
389 const Double_t *pid=track->PID();
392 for (Int_t i=0; i<10; i++) {
393 if (pid[i]>max) {k=i; max=pid[i];}
396 if(k==2) type[0]=kTRUE;
397 if(k==3) type[1]=kTRUE;
398 if(k==4) type[2]=kTRUE;
402 //--------------------------------
403 Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
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;
409 //--------------------------------
410 Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
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;
418 //--------------------------------
419 Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
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;
429 //--------------------------------
430 Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
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;
436 //--------------------------------
437 Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
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);
444 AliError("Wrong detector name");
448 //--------------------------------------------
449 Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
450 // TPC nsigma cut PID, different sigmas in different p bins
452 AliAODPid *pidObj = track->GetDetPid();
453 Double_t mom = pidObj->GetTPCmomentum();
454 if(mom>fPtThresholdTPC) return kTRUE;
457 if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
458 nsigma=TMath::Abs(nsigma);
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;
468 Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
469 // combination of the PID info coming from TPC and TOF
471 Double_t ptrack=track->P();
472 if(ptrack>fMaxTrackMomForCombinedPID) return 1;
474 Bool_t okTPC=CheckTPCPIDStatus(track);
475 if(ptrack>fPtThresholdTPC) okTPC=kFALSE;
476 Bool_t okTOF=CheckTOFPIDStatus(track);
479 //TOF || TPC (a la' Andrea R.)
481 // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
482 // the method returns the sum of the response of the 2 detectors
485 if(!okTPC && !okTOF) return 0;
492 if(TPCRawAsym(track,specie)) tTPCinfo=1;
494 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
496 if(fCompat && tTPCinfo<0){
497 Double_t sig0tmp=fnSigma[0];
498 SetSigma(0,fnSigmaCompat[0]);
499 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
506 if(!okTOF && fTPC) return tTPCinfo;
508 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
509 if(fCompat && tTOFinfo>0){
510 if(ptrack>fPCompatTOF) {
511 if(ApplyTOFCompatibilityBand(track,specie)==specie) tTOFinfo=0;
517 if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
518 if(!okTOF) return tTPCinfo;
522 if(tTPCinfo+tTOFinfo==0 && fITS){
523 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
525 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
528 return tTPCinfo+tTOFinfo;
532 //TPC & TOF (a la' Yifei)
533 // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
539 if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
541 if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
547 if(fTPC && !okTOF) return tTPCinfo;
548 if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
551 if(tTOFinfo==1 && tTPCinfo==1) return 1;
553 if(tTPCinfo+tTOFinfo==0 && fITS){
554 if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
556 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
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;
569 if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
572 if(TPCRawAsym(track,specie)) tTPCinfo=1;
574 if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
580 if(ptrack>=fPLimit[1] && fTOF){
582 if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
587 if(ptrack<fPLimit[0] && fITS){
588 if(!CheckITSPIDStatus(track)) return 0;
589 if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
594 if(fMatch==4 || fMatch==5){
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
601 Double_t nSigmaTPC=0.;
603 nSigmaTPC=fPidResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)specie);
604 if(nSigmaTPC<-990.) nSigmaTPC=0.;
606 Double_t nSigmaTOF=0.;
608 nSigmaTOF=fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)specie);
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;
613 Double_t nSigma2=nSigmaTPC*nSigmaTPC+nSigmaTOF*nSigmaTOF;
614 if(nSigma2<fMaxnSigmaCombined[iPart]*fMaxnSigmaCombined[iPart]) return 1;
618 if((nSigmaTPC>fMinnSigmaTPC[iPart] && nSigmaTPC<fMaxnSigmaTPC[iPart]) &&
619 (nSigmaTOF>fMinnSigmaTOF[iPart] && nSigmaTOF<fMaxnSigmaTOF[iPart])) return 1;
628 //----------------------------------
629 Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
630 // general method to compute PID
632 return MatchTPCTOF(track,specie);
634 if(fTPC && !fTOF && !fITS) {
637 tTPCres=ApplyPidTPCRaw(track,specie);
638 if(tTPCres==specie) return 1;
641 if(TPCRawAsym(track,specie)) tTPCres=1;
645 }else if(fTOF && !fTPC && !fITS) {
646 Int_t tTOFres=ApplyPidTOFRaw(track,specie);
647 if(tTOFres==specie) return 1;
649 }else if(fITS && !fTPC && !fTOF) {
650 Int_t tITSres=ApplyPidITSRaw(track,specie);
651 if(tITSres==specie) return 1;
654 AliError("You should enable just one detector if you don't want to match");
659 //--------------------------------------------
660 void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
661 // TPC bethe bloch parameters
664 if(fPbPb) { // PbPb MC
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;
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;
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;
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;
694 } else { // Real Data
696 if(fOnePad) { // pp 1-pad (since LHC10d)
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;
704 } else if(fPbPb) { // PbPb
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;
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;
718 } else if(fppLowEn2011){ // pp low energy
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;
726 } else { // pp no 1-pad (LHC10bc)
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;
740 //-----------------------
741 void AliAODPidHF::SetBetheBloch() {
742 // Set Bethe Bloch Parameters
744 Double_t alephParameters[5];
745 GetTPCBetheBlochParams(alephParameters);
746 fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
752 //--------------------------------------------------------------------------
753 Int_t AliAODPidHF::GetnSigmaITS(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
754 // get n sigma for ITS
757 if (!CheckITSPIDStatus(track)) return -1;
759 Double_t nsigmaITS=-999;
762 Double_t mom=track->P();
763 AliAODPid *pidObj = track->GetDetPid();
764 Double_t dedx=pidObj->GetITSsignal();
766 AliITSPIDResponse itsResponse;
767 AliPID::EParticleType type=AliPID::EParticleType(species);
768 nsigmaITS = itsResponse.GetNumberOfSigmas(mom,dedx,type);
773 AliPID::EParticleType type=AliPID::EParticleType(species);
774 nsigmaITS = fPidResponse->NumberOfSigmasITS(track,type);
783 //--------------------------------------------------------------------------
784 Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
785 // get n sigma for TPC
787 if(!CheckTPCPIDStatus(track)) return -1;
789 Double_t nsigmaTPC=-999;
792 AliAODPid *pidObj = track->GetDetPid();
793 Double_t dedx=pidObj->GetTPCsignal();
794 Double_t mom = pidObj->GetTPCmomentum();
795 if(mom>fPtThresholdTPC) return -2;
796 UShort_t nTPCClus=pidObj->GetTPCsignalN();
797 if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
798 AliPID::EParticleType type=AliPID::EParticleType(species);
799 nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
802 if(!fPidResponse) return -1;
803 AliPID::EParticleType type=AliPID::EParticleType(species);
804 nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
810 //-----------------------------
812 Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
813 // get n sigma for TOF
815 if(!CheckTOFPIDStatus(track)) return -1;
818 nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
821 AliFatal("To use TOF PID you need to attach AliPIDResponseTask");
827 //-----------------------
828 Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
829 // Exclude a given hypothesis (labelTracks) in detector
831 if (detectors.Contains("ITS")) {
833 AliInfo("Nothing to be done");
836 if (GetnSigmaITS(track,labelTrack,nsigma)==1){
837 if(nsigma>nsigmaCut) return kTRUE;
842 } else if (detectors.Contains("TPC")) {
845 if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
846 if(nsigma>nsigmaCut) return kTRUE;
850 } else if (detectors.Contains("TOF")) {
853 if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
854 if(nsigma>nsigmaCut) return kTRUE;
862 //-----------------------
863 Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
864 // TOF proton compatibility
866 if(!CheckTOFPIDStatus(track)) return 0;
869 if(GetnSigmaTOF(track,3,nsigma)==1){
870 if(nsigma>nsigmaK) return kTRUE;
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();
879 AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
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));
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
893 if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
899 //--------------------------------------------------------------------------
900 void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
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
907 GetPidCombined()->SetPriorDistribution(type,prior);
909 //--------------------------------------------------------------------------
910 void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
913 // Drawing prior distribution for type "type"
916 GetPidCombined()->GetPriorDistribution(type)->Draw();
919 //-----------------------------
920 void AliAODPidHF::SetPriorsHistos(TString priorFileName){
921 // Set histograms with priors
923 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
924 if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
927 nt+=AliPID::ParticleName(ispecies);
929 TDirectory *current = gDirectory;
930 TFile *priorFile=TFile::Open(priorFileName);
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"));
936 fPriorsH[AliPID::kProton] = new TH1F(*h3);
937 fPriorsH[AliPID::kKaon ] = new TH1F(*h2);
938 fPriorsH[AliPID::kPion ] = new TH1F(*h1);
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);
948 //----------------------------------
949 void AliAODPidHF::SetUpCombinedPID(){
950 // Configuration of combined Bayesian PID
952 fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
954 for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
955 fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
958 fPidCombined->SetDefaultTPCPriors();
960 switch (fCombDetectors){
962 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
965 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC & AliPIDResponse::kDetTOF);
968 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
971 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
974 fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
980 //-----------------------------
981 void 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 ");
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);
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);
1001 printf("--- Matching algorithm = %d ---\n",fMatch);
1003 if(fITS) printf("nSigmaITS = %.2f\n",fnSigma[4]);
1005 printf("nSigmaTOF = %.2f\n",fnSigma[3]);
1006 if(fCompat) printf("Compatibility band at nSigmaTOF=%.2f for p>%.2f\n",fnSigmaCompat[1],fPCompatTOF);
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]);
1015 printf("nSigmaTPC = %.2f\n",fnSigma[0]);
1017 if(fCompat) printf("Compatibility band at nSigmaTPC=%.2f\n",fnSigmaCompat[0]);
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]);