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