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(kFALSE),
  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   Double_t nsigma;
 217   Int_t pid=-1;
 218
 219   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
 220     Double_t nsigmaMin=999.;
 221     for(Int_t ipart=0;ipart<5;ipart++){
 222       if(GetnSigmaTPC(track,ipart,nsigma)==1){
 223         nsigma=TMath::Abs(nsigma);
 224         if((nsigma<nsigmaMin) && (nsigma<fnSigma[0])) {
 225           pid=ipart;
 226           nsigmaMin=nsigma;
 227         }
 228       }
 229     }
 230   }else{ // asks only for one particle specie
 231     if(GetnSigmaTPC(track,specie,nsigma)==1){
 232       nsigma=TMath::Abs(nsigma);
 233       if (nsigma>fnSigma[0]) pid=-1;
 234       else pid=specie;
 235     }
 236   }
 237
 238   return pid;
 239 }
 240 //----------------------------
 241 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
 242 // truncated mean, ITS PID
 243
 244   if(!CheckITSPIDStatus(track)) return 0;
 245   Int_t pid=-1;
 246
 247   if(fOldPid){
 248   Double_t mom=track->P();
 249   AliAODPid *pidObj = track->GetDetPid();
 250
 251   Double_t dedx=pidObj->GetITSsignal();
 252   UChar_t clumap=track->GetITSClusterMap();
 253   Int_t nPointsForPid=0;
 254   for(Int_t i=2; i<6; i++){
 255    if(clumap&(1<<i)) ++nPointsForPid;
 256   }
 257
 258   Bool_t isSA=kTRUE;
 259   if(track->GetStatus() & AliESDtrack::kTPCin) isSA = kFALSE;
 260
 261   AliITSPIDResponse itsResponse;
 262   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
 263    Double_t nsigmaMax=fnSigma[4];
 264    for(Int_t ipart=0;ipart<5;ipart++){
 265     AliPID::EParticleType type=AliPID::EParticleType(ipart);
 266     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type,nPointsForPid,isSA));
 267     if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
 268      pid=ipart;
 269      nsigmaMax=nsigma;
 270     }
 271    }
 272   }else{ // asks only for one particle specie
 273    AliPID::EParticleType type=AliPID::EParticleType(specie);
 274     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
 275    if (nsigma>fnSigma[4]) {
 276     pid=-1;
 277    }else{
 278     pid=specie;
 279    }
 280   }
 281  }else{ // old pid
 282
 283   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
 284    Double_t nsigmaMax=fnSigma[4];
 285    for(Int_t ipart=0;ipart<5;ipart++){
 286     AliPID::EParticleType type=AliPID::EParticleType(ipart);
 287     Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasITS(track,type));
 288     if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
 289      pid=ipart;
 290      nsigmaMax=nsigma;
 291     }
 292    }
 293   }else{ // asks only for one particle specie
 294    AliPID::EParticleType type=AliPID::EParticleType(specie);
 295     Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasITS(track,type));
 296    if (nsigma>fnSigma[4]) {
 297     pid=-1;
 298    }else{
 299     pid=specie;
 300    }
 301   }
 302  } //new pid
 303
 304  return pid;
 305 }
 306 //----------------------------
 307 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
 308 // n-sigma cut, TOF PID
 309
 310   Double_t nsigma;
 311   Int_t pid=-1;
 312
 313   if(specie<0){
 314     Double_t nsigmaMin=999.;
 315     for(Int_t ipart=0;ipart<5;ipart++){
 316       if(GetnSigmaTOF(track,ipart,nsigma)==1){
 317         nsigma=TMath::Abs(nsigma);
 318         if((nsigma<nsigmaMin)&& (nsigma<fnSigma[3])){
 319           pid=ipart;
 320           nsigmaMin=nsigma;
 321         }
 322       }
 323     }
 324   }else{ // asks only for one particle specie
 325     if(GetnSigmaTOF(track,specie,nsigma)==1){
 326       nsigma=TMath::Abs(nsigma);
 327       if (nsigma>fnSigma[3]) pid=-1;
 328       else pid=specie;
 329     }
 330   }
 331   return pid;
 332   /*
 333  Double_t time[AliPID::kSPECIESN];
 334  Double_t sigmaTOFPid[AliPID::kSPECIES];
 335  AliAODPid *pidObj = track->GetDetPid();
 336  pidObj->GetIntegratedTimes(time);
 337  Double_t sigTOF=pidObj->GetTOFsignal();
 338
 339  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
 340  if (event) {
 341    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
 342    if (tofH && fPidResponse) { // reading new AOD with new aliroot
 343      AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
 344      sigTOF -= TOFres.GetStartTime(track->P());
 345      if (specie<0) {
 346        for (Int_t ipart = 0; ipart<5; ipart++) {
 347          sigmaTOFPid[ipart]=TOFres.GetExpectedSigma(track->P(),time[ipart],AliPID::ParticleMass(ipart));
 348        }
 349      }
 350      else sigmaTOFPid[specie]=TOFres.GetExpectedSigma(track->P(),time[specie],AliPID::ParticleMass(specie)); //fTOFResponse is set in InitialiseEvent
 351    } else  pidObj->GetTOFpidResolution(sigmaTOFPid); // reading old AOD with new aliroot
 352  } else  pidObj->GetTOFpidResolution(sigmaTOFPid);  //reading old AOD with old aliroot
 353
 354  Int_t pid=-1;
 355
 356   if(specie<0){
 357    Double_t sigmaTOFtrack;
 358    if (sigmaTOFPid[4]>0) sigmaTOFtrack=sigmaTOFPid[4];
 359    else sigmaTOFtrack=fTOFSigma;
 360    Double_t nsigmaMax=sigmaTOFtrack*fnSigma[3];
 361    for(Int_t ipart=0;ipart<5;ipart++){
 362     Double_t nsigma=TMath::Abs(sigTOF-time[ipart]);
 363     if (sigmaTOFPid[ipart]>0) sigmaTOFtrack=sigmaTOFPid[ipart];
 364     else sigmaTOFtrack=fTOFSigma;  // backward compatibility for old AODs
 365     if((nsigma<nsigmaMax) && (nsigma<fnSigma[3]*sigmaTOFtrack)) {
 366      pid=ipart;
 367      nsigmaMax=nsigma;
 368     }
 369    }
 370   }else{ // asks only for one particle specie
 371     Double_t nsigma=TMath::Abs(sigTOF-time[specie]);
 372     Double_t sigmaTOFtrack;
 373     if (sigmaTOFPid[specie]>0) sigmaTOFtrack=sigmaTOFPid[specie];
 374     else sigmaTOFtrack=fTOFSigma;  // backward compatibility for old AODs
 375     if (nsigma>fnSigma[3]*sigmaTOFtrack) {
 376       pid=-1;
 377     }else{
 378       pid=specie;
 379     }
 380   }
 381  return pid;
 382   */
 383 }
 384 //------------------------------
 385 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
 386 // combined PID stored inside the AOD track
 387
 388  const Double_t *pid=track->PID();
 389  Float_t max=0.;
 390  Int_t k=-1;
 391  for (Int_t i=0; i<10; i++) {
 392    if (pid[i]>max) {k=i; max=pid[i];}
 393  }
 394
 395  if(k==2) type[0]=kTRUE;
 396  if(k==3) type[1]=kTRUE;
 397  if(k==4) type[2]=kTRUE;
 398
 399  return;
 400 }
 401 //--------------------
 402 void AliAODPidHF::BayesianProbability(AliAODTrack *track,Double_t *pid) const{
 403 // bayesian PID for single detectors or combined
 404
 405   if(fITS && !fTPC && !fTOF) {BayesianProbabilityITS(track,pid);return;}
 406   if(fTPC && !fITS && !fTOF) {BayesianProbabilityTPC(track,pid);return;}
 407   if(fTOF && !fITS && !fTPC) {BayesianProbabilityTOF(track,pid);return;}
 408
 409     Double_t probITS[5]={1.,1.,1.,1.,1.};
 410     Double_t probTPC[5]={1.,1.,1.,1.,1.};
 411     Double_t probTOF[5]={1.,1.,1.,1.,1.};
 412     if(fITS) BayesianProbabilityITS(track,probITS);
 413     if(fTPC) BayesianProbabilityTPC(track,probTPC);
 414     if(fTOF) BayesianProbabilityTOF(track,probTOF);
 415     Double_t probTot[5]={0.,0.,0.,0.,0.};
 416     for(Int_t i=0;i<5;i++){
 417      probTot[i]=probITS[i]*probTPC[i]*probTOF[i];
 418     }
 419     for(Int_t i2=0;i2<5;i2++){
 420      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]);
 421     }
 422
 423  return;
 424
 425 }
 426 //------------------------------------
 427 void AliAODPidHF::BayesianProbabilityITS(AliAODTrack *track,Double_t *prob) const{
 428
 429 // bayesian PID for ITS
 430  AliAODpidUtil pid;
 431  Double_t itspid[AliPID::kSPECIES];
 432  pid.MakeITSPID(track,itspid);
 433  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 434   if(fTOF || fTPC || fTRD){
 435    prob[ind]=itspid[ind];
 436   }else{
 437    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]);
 438   }
 439  }
 440  return;
 441
 442 }
 443 //------------------------------------
 444 void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack *track,Double_t *prob) const{
 445 // bayesian PID for TPC
 446
 447  AliAODpidUtil pid;
 448  Double_t tpcpid[AliPID::kSPECIES];
 449  pid.MakeTPCPID(track,tpcpid);
 450  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 451   if(fTOF || fITS || fTRD){
 452    prob[ind]=tpcpid[ind];
 453   }else{
 454    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]);
 455  }
 456 }
 457  return;
 458
 459 }
 460 //------------------------------------
 461 void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack *track,Double_t *prob) const{
 462 // bayesian PID for TOF
 463
 464  AliAODpidUtil pid;
 465  Double_t tofpid[AliPID::kSPECIES];
 466  pid.MakeTOFPID(track,tofpid);
 467  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 468   if(fTPC || fITS || fTRD){
 469    prob[ind]=tofpid[ind];
 470   }else{
 471   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]);
 472  }
 473 }
 474  return;
 475
 476 }
 477 //---------------------------------
 478 void AliAODPidHF::BayesianProbabilityTRD(AliAODTrack *track,Double_t *prob) const{
 479 // bayesian PID for TRD
 480
 481  AliAODpidUtil pid;
 482  Double_t trdpid[AliPID::kSPECIES];
 483  pid.MakeTRDPID(track,trdpid);
 484  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 485   if(fTPC || fITS || fTOF){
 486    prob[ind]=trdpid[ind];
 487   }else{
 488    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]);
 489  }
 490 }
 491   return;
 492
 493  }
 494 //--------------------------------
 495 Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
 496   // ITS PID quality cuts
 497   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
 498   UChar_t clumap=track->GetITSClusterMap();
 499   Int_t nPointsForPid=0;
 500   for(Int_t i=2; i<6; i++){
 501     if(clumap&(1<<i)) ++nPointsForPid;
 502   }
 503   if(nPointsForPid<3) return kFALSE;
 504   return kTRUE;
 505 }
 506 //--------------------------------
 507 Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
 508   // TPC PID quality cuts
 509   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
 510   UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
 511   if (nTPCClus<70) return kFALSE;
 512   return kTRUE;
 513 }
 514 //--------------------------------
 515 Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
 516   // TOF PID quality cuts
 517   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return kFALSE;
 518   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return kFALSE;
 519   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return kFALSE;
 520   if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0)    return kFALSE;
 521   return kTRUE;
 522 }
 523 //--------------------------------
 524 Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
 525   // TRD PID quality cuts
 526   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return kFALSE;
 527   return kTRUE;
 528 }
 529 //--------------------------------
 530 Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
 531
 532 // Quality cuts on the tracks, detector by detector
 533
 534  if(detectors.Contains("ITS")){
 535   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
 536   UChar_t clumap=track->GetITSClusterMap();
 537   Int_t nPointsForPid=0;
 538   for(Int_t i=2; i<6; i++){
 539    if(clumap&(1<<i)) ++nPointsForPid;
 540   }
 541   if(nPointsForPid<3) return kFALSE;
 542  }
 543
 544  if(detectors.Contains("TPC")){
 545    if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
 546    UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
 547    if (nTPCClus<70) return kFALSE;
 548  }
 549
 550  if(detectors.Contains("TOF")){
 551    if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return kFALSE;
 552    if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return kFALSE;
 553    if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return kFALSE;
 554    if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0)    return kFALSE;
 555  }
 556
 557
 558  if(detectors.Contains("TRD")){
 559   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return kFALSE;
 560   //UChar_t ntracklets = track->GetTRDntrackletsPID();
 561   //if(ntracklets<4) return kFALSE;
 562  }
 563
 564  return kTRUE;
 565 }
 566 //--------------------------------------------
 567 Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
 568 // TPC nsigma cut PID, different sigmas in different p bins
 569
 570   Double_t nsigma;
 571   if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
 572   nsigma=TMath::Abs(nsigma);
 573
 574   AliAODPid *pidObj = track->GetDetPid();
 575   Double_t mom = pidObj->GetTPCmomentum();
 576   if(mom>fPtThresholdTPC) return 0;
 577
 578   if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
 579   if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
 580   if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
 581
 582   return kFALSE;
 583 }
 584 //------------------
 585 Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
 586   // combination of the PID info coming from TPC and TOF
 587
 588   Bool_t okTPC=CheckTPCPIDStatus(track);
 589   Bool_t okTOF=CheckTOFPIDStatus(track);
 590
 591   if(fMatch==1){
 592     //TOF || TPC (a la' Andrea R.)
 593     // convention:
 594     // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
 595     // the method returns the sum of the response of the 2 detectors
 596
 597     if(fTPC && fTOF) {
 598       if(!okTPC && !okTOF) return 0;
 599     }
 600
 601     Int_t tTPCinfo=0;
 602     if(fTPC && okTPC){
 603       tTPCinfo=-1;
 604       if(fAsym) {
 605         if(TPCRawAsym(track,specie)) tTPCinfo=1;
 606       }else{
 607         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
 608       }
 609       if(fCompat && tTPCinfo<0){
 610         Double_t sig0tmp=fnSigma[0];
 611         SetSigma(0,fnSigmaCompat[0]);
 612         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
 613         SetSigma(0,sig0tmp);
 614       }
 615     }
 616
 617     Int_t tTOFinfo=0;
 618     if(fTOF){
 619       if(!okTOF && fTPC) return tTPCinfo;
 620       tTOFinfo=-1;
 621       if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
 622       if(fCompat && tTOFinfo>0){
 623         Double_t ptrack=track->P();
 624         if(ptrack>fPCompatTOF) {
 625           Double_t sig0tmp=fnSigma[3];
 626           SetSigma(3,fnSigmaCompat[1]);
 627           if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=0;
 628           SetSigma(3,sig0tmp);
 629         }
 630       }
 631     }
 632
 633
 634     if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
 635       if(!okTOF) return tTPCinfo;
 636       return tTOFinfo;
 637     }
 638
 639     if(tTPCinfo+tTOFinfo==0 && fITS){
 640       if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
 641       Int_t tITSinfo = -1;
 642       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 643       return tITSinfo;
 644     }
 645     return tTPCinfo+tTOFinfo;
 646   }
 647
 648   if(fMatch==2){
 649     //TPC & TOF (a la' Yifei)
 650     // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
 651     Int_t tTPCinfo=0;
 652
 653     if(fTPC && okTPC) {
 654       tTPCinfo=1;
 655       if(fAsym){
 656         if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
 657       }else{
 658         if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
 659       }
 660     }
 661
 662     Int_t tTOFinfo=1;
 663     if(fTOF){
 664       if(fTPC && !okTOF) return tTPCinfo;
 665       if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
 666     }
 667
 668     if(tTOFinfo==1 && tTPCinfo==1) return 1;
 669
 670     if(tTPCinfo+tTOFinfo==0 && fITS){
 671       if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
 672       Int_t tITSinfo = -1;
 673       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 674       return tITSinfo;
 675     }
 676     return -1;
 677   }
 678
 679   if(fMatch==3){
 680     //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
 681     // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
 682     if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
 683
 684     Double_t ptrack=track->P();
 685
 686     Int_t tTPCinfo=-1;
 687     if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
 688       if(!okTPC) return 0;
 689       if(fAsym) {
 690         if(TPCRawAsym(track,specie)) tTPCinfo=1;
 691       }else{
 692         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
 693       }
 694       return tTPCinfo;
 695     }
 696
 697     Int_t tTOFinfo=-1;
 698     if(ptrack>=fPLimit[1] && fTOF){
 699       if(!okTOF) return 0;
 700       if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
 701       return tTOFinfo;
 702     }
 703
 704     Int_t tITSinfo=-1;
 705     if(ptrack<fPLimit[0] && fITS){
 706       if(!CheckITSPIDStatus(track)) return 0;
 707       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 708       return tITSinfo;
 709     }
 710   }
 711
 712   return -1;
 713
 714 }
 715 //----------------------------------
 716 Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
 717 // general method to compute PID
 718   if(fMatch>0){
 719     return MatchTPCTOF(track,specie);
 720   }else{
 721     if(fTPC && !fTOF && !fITS) {
 722       Int_t tTPCres=0;
 723       if(!fAsym){
 724         tTPCres=ApplyPidTPCRaw(track,specie);
 725         if(tTPCres==specie) return 1;
 726         else return tTPCres;
 727       }else{
 728         if(TPCRawAsym(track,specie)) tTPCres=1;
 729         else tTPCres=-1;
 730       }
 731       return tTPCres;
 732     }else if(fTOF && !fTPC && !fITS) {
 733       Int_t tTOFres=ApplyPidTOFRaw(track,specie);
 734       if(tTOFres==specie) return 1;
 735       else return tTOFres;
 736     }else if(fITS && !fTPC && !fTOF) {
 737       Int_t tITSres=ApplyPidITSRaw(track,specie);
 738       if(tITSres==specie) return 1;
 739       else return tITSres;
 740     }else{
 741       AliError("You should enable just one detector if you don't want to match");
 742       return 0;
 743     }
 744   }
 745 }
 746 //--------------------------------------------
 747 void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
 748   // TPC bethe bloch parameters
 749  if(fMC) {  // MC
 750
 751    if(fPbPb) { // PbPb MC
 752
 753      alephParameters[0] = 1.44405/50.;
 754      alephParameters[1] = 2.35409e+01;
 755      alephParameters[2] = TMath::Exp(-2.90330e+01);
 756      alephParameters[3] = 2.10681e+00;
 757      alephParameters[4] = 4.62254e+00;
 758
 759    } else {  // pp MC
 760      if(fMCLowEn2011){
 761        alephParameters[0]=0.0207667;
 762        alephParameters[1]=29.9936;
 763        alephParameters[2]=3.87866e-11;
 764        alephParameters[3]=2.17291;
 765        alephParameters[4]=7.1623;
 766      }else if(fOnePad){
 767        alephParameters[0]=0.029021;
 768        alephParameters[1]=25.4181;
 769        alephParameters[2]=4.66596e-08;
 770        alephParameters[3]=1.90008;
 771        alephParameters[4]=4.63783;
 772      }else{
 773        alephParameters[0] = 2.15898/50.;
 774        alephParameters[1] = 1.75295e+01;
 775        alephParameters[2] = 3.40030e-09;
 776        alephParameters[3] = 1.96178e+00;
 777        alephParameters[4] = 3.91720e+00;
 778      }
 779    }
 780
 781  } else { // Real Data
 782
 783    if(fOnePad) { // pp 1-pad (since LHC10d)
 784
 785      alephParameters[0] =1.34490e+00/50.;
 786      alephParameters[1] = 2.69455e+01;
 787      alephParameters[2] = TMath::Exp(-2.97552e+01);
 788      alephParameters[3] = 2.35339e+00;
 789      alephParameters[4] = 5.98079e+00;
 790
 791    } else if(fPbPb) { // PbPb
 792
 793      // alephParameters[0] = 1.25202/50.;
 794      // alephParameters[1] = 2.74992e+01;
 795      // alephParameters[2] = TMath::Exp(-3.31517e+01);
 796      // alephParameters[3] = 2.46246;
 797      // alephParameters[4] = 6.78938;
 798
 799      alephParameters[0] = 5.10207e+00/50.;
 800      alephParameters[1] = 7.94982e+00;
 801      alephParameters[2] = TMath::Exp(-9.07942e+00);
 802      alephParameters[3] = 2.38808e+00;
 803      alephParameters[4] = 1.68165e+00;
 804
 805    } else if(fppLowEn2011){ // pp low energy
 806
 807      alephParameters[0]=0.031642;
 808      alephParameters[1]=22.353;
 809      alephParameters[2]=4.16239e-12;
 810      alephParameters[3]=2.61952;
 811      alephParameters[4]=5.76086;
 812
 813    } else {  // pp no 1-pad (LHC10bc)
 814
 815      alephParameters[0] = 0.0283086/0.97;
 816      alephParameters[1] = 2.63394e+01;
 817      alephParameters[2] = 5.04114e-11;
 818      alephParameters[3] = 2.12543e+00;
 819      alephParameters[4] = 4.88663e+00;
 820
 821    }
 822
 823  }
 824
 825 }
 826
 827 //-----------------------
 828 void AliAODPidHF::SetBetheBloch() {
 829   // Set Bethe Bloch Parameters
 830
 831  Double_t alephParameters[5];
 832  GetTPCBetheBlochParams(alephParameters);
 833  fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 834
 835  return;
 836 }
 837 //-----------------------
 838 Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
 839
 840   if(!CheckTOFPIDStatus(track)) return 0;
 841
 842   Double_t nsigma;
 843   if(GetnSigmaTOF(track,3,nsigma)==1){
 844     if(nsigma>nsigmaK) return kTRUE;
 845   }
 846   return kFALSE;
 847   /*  Double_t time[AliPID::kSPECIESN];
 848   Double_t sigmaTOFPid[AliPID::kSPECIES];
 849   AliAODPid *pidObj = track->GetDetPid();
 850   pidObj->GetIntegratedTimes(time);
 851   Double_t sigTOF=pidObj->GetTOFsignal();
 852
 853  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
 854  if (event) {
 855    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
 856    if (tofH && fPidResponse) {
 857      AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
 858      sigTOF -= TOFres.GetStartTime(track->P());
 859      sigmaTOFPid[3]=TOFres.GetExpectedSigma(track->P(),time[3],AliPID::ParticleMass(3));
 860    }
 861    else  pidObj->GetTOFpidResolution(sigmaTOFPid);
 862  } else  pidObj->GetTOFpidResolution(sigmaTOFPid);
 863   Double_t sigmaTOFtrack;
 864   if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
 865   else sigmaTOFtrack=fTOFSigma;  // backward compatibility for old AODs
 866
 867   if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
 868
 869   return kFALSE;
 870   */
 871 }
 872 //--------------------------------------------------------------------------
 873 void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
 874
 875         //
 876         // method setting the prior distributions to the AliPIDCombined object of the AliAODPidHF data member
 877         // all the checks are done directly in the AliPIDCombined object
 878         //
 879
 880         GetPidCombined()->SetPriorDistribution(type,prior);
 881 }
 882 //--------------------------------------------------------------------------
 883 void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
 884
 885         //
 886         // Drawing prior distribution for type "type"
 887
 888         new TCanvas();
 889         GetPidCombined()->GetPriorDistribution(type)->Draw();
 890 }
 891
 892 //--------------------------------------------------------------------------
 893 Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
 894
 895   if(!CheckTPCPIDStatus(track)) return -1;
 896
 897   Double_t nsigmaTPC=-999;
 898
 899   if(fOldPid){
 900     AliAODPid *pidObj = track->GetDetPid();
 901     Double_t dedx=pidObj->GetTPCsignal();
 902     Double_t mom = pidObj->GetTPCmomentum();
 903     if(mom>fPtThresholdTPC) return -2;
 904     UShort_t nTPCClus=pidObj->GetTPCsignalN();
 905     if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
 906     AliPID::EParticleType type=AliPID::EParticleType(species);
 907     nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
 908     nsigma=nsigmaTPC;
 909   } else{
 910     if(!fPidResponse) return -1;
 911     AliPID::EParticleType type=AliPID::EParticleType(species);
 912     nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
 913     nsigma=nsigmaTPC;
 914   }
 915   return 1;
 916 }
 917
 918 //-----------------------------
 919
 920 Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
 921
 922   if(!CheckTOFPIDStatus(track)) return -1;
 923
 924   if(fPidResponse){
 925     nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
 926   }else{
 927     AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
 928     if (event) {
 929       AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
 930       if (tofH) { // reading a new AOD without setting fPidResponse!!!
 931         AliFatal("To use this AOD you must start AliPIDResponseTask");
 932       }
 933       else {
 934         Double_t time[AliPID::kSPECIESN];
 935         Double_t sigmaTOFPid[AliPID::kSPECIES];
 936         AliAODPid *pidObj = track->GetDetPid();
 937         pidObj->GetIntegratedTimes(time);
 938         Double_t sigTOF=pidObj->GetTOFsignal();
 939         pidObj->GetTOFpidResolution(sigmaTOFPid);
 940         if(sigmaTOFPid[species]<1e-99) return -2;
 941         Double_t sigmaTOF=(sigTOF-time[species])/sigmaTOFPid[species];
 942         nsigma=sigmaTOF;
 943       }
 944     }
 945   }
 946   return 1;
 947 }
 948
 949 //-----------------------
 950 Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut,
 951                                TString detectors) {
 952
 953   if (detectors.Contains("ITS")) {
 954
 955     AliInfo("Nothing to be done");
 956     /*
 957     Double_t nsigma=0.;
 958     if (GetnSigmaITS(track,labelTrack,nsigma)==1){
 959       if(nsigma>nsigmaCut) return kTRUE;
 960     }
 961     */
 962     return kFALSE;
 963
 964   } else if (detectors.Contains("TPC")) {
 965
 966     Double_t nsigma=0.;
 967     if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
 968       if(nsigma>nsigmaCut) return kTRUE;
 969     }
 970     return kFALSE;
 971
 972   } else if (detectors.Contains("TOF")) {
 973
 974     if (!(CheckTOFPIDStatus(track))) return kFALSE;
 975     Double_t nsigma=0.;
 976     if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
 977       if(nsigma>nsigmaCut) return kTRUE;
 978     }
 979     return kFALSE;
 980
 981   }
 982   return kFALSE;
 983
 984 }
 985 //-----------------------------