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