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