PWG3/vertexingHF/AliAODPidHF.cxx

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