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