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