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