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, J. van der Maarel j.vandermaarel@cern.ch
  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 fForceTOFforKaons(kFALSE),
  59 fCompat(kFALSE),
  60 fPCompatTOF(1.5),
  61 fUseAsymTOF(kFALSE),
  62 fLownSigmaTOF(-3.),
  63 fUpnSigmaTOF(3.),
  64 fLownSigmaCompatTOF(-3.),
  65 fUpnSigmaCompatTOF(3.),
  66 fnNSigmaCompat(2),
  67 fnSigmaCompat(0),
  68 fMC(kFALSE),
  69 fOnePad(kFALSE),
  70 fMCLowEn2011(kFALSE),
  71 fppLowEn2011(kFALSE),
  72 fPbPb(kFALSE),
  73 fTOFdecide(kFALSE),
  74 fOldPid(kFALSE),
  75 fPtThresholdTPC(999999.),
  76 fMaxTrackMomForCombinedPID(999999.),
  77 fPidResponse(0),
  78 fPidCombined(new AliPIDCombined()),
  79 fTPCResponse(new AliTPCPIDResponse()),
  80 fPriorsH(),
  81 fCombDetectors(kTPCTOF),
  82 fUseCombined(kFALSE),
  83 fDefaultPriors(kTRUE)
  84 {
  85   //
  86   // Default constructor
  87   //
  88   fPLimit=new Double_t[fnPLimit];
  89   fnSigma=new Double_t[fnNSigma];
  90   fPriors=new Double_t[fnPriors];
  91   fnSigmaCompat=new Double_t[fnNSigmaCompat];
  92
  93   for(Int_t i=0;i<fnNSigma;i++){
  94     fnSigma[i]=0.;
  95   }
  96   for(Int_t i=0;i<fnPriors;i++){
  97     fPriors[i]=0.;
  98   }
  99   for(Int_t i=0;i<fnPLimit;i++){
 100     fPLimit[i]=0.;
 101   }
 102   for(Int_t i=0;i<fnNSigmaCompat;i++){
 103     fnSigmaCompat[i]=3.;
 104   }
 105   for(Int_t i=0; i<3; i++){ // pi, K, proton
 106     fMaxnSigmaCombined[i]=3.;
 107     fMinnSigmaTPC[i]=-3;
 108     fMaxnSigmaTPC[i]=3;
 109     fMinnSigmaTOF[i]=-3;
 110     fMaxnSigmaTOF[i]=3;
 111   }
 112   for (Int_t s=0;s<AliPID::kSPECIES;s++) {
 113     for (Int_t d=0;d<4;d++) {
 114       fIdBandMin[s][d] = NULL;
 115       fIdBandMax[s][d] = NULL;
 116       fCompBandMin[s][d] = NULL;
 117       fCompBandMax[s][d] = NULL;
 118     }
 119   }
 120
 121 }
 122 //----------------------
 123 AliAODPidHF::~AliAODPidHF()
 124 {
 125   // destructor
 126   if(fPLimit) delete [] fPLimit;
 127   if(fnSigma) delete [] fnSigma;
 128   if(fPriors) delete [] fPriors;
 129   if(fnSigmaCompat) delete [] fnSigmaCompat;
 130   delete fPidCombined;
 131
 132   delete fTPCResponse;
 133   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
 134     delete fPriorsH[ispecies];
 135   }
 136
 137   for (Int_t s=0;s<AliPID::kSPECIES;s++) {
 138     for (Int_t d=0;d<4;d++) {
 139       delete fIdBandMin[s][d];
 140       delete fIdBandMax[s][d];
 141       delete fCompBandMin[s][d];
 142       delete fCompBandMax[s][d];
 143     }
 144   }
 145 }
 146 //------------------------
 147 AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
 148 TObject(),
 149 fnNSigma(pid.fnNSigma),
 150 fnSigma(0),
 151 fTOFSigma(pid.fTOFSigma),
 152 fCutTOFmismatch(pid.fCutTOFmismatch),
 153 fMinNClustersTPCPID(pid.fMinNClustersTPCPID),
 154 fnPriors(pid.fnPriors),
 155 fPriors(0),
 156 fnPLimit(pid.fnPLimit),
 157 fPLimit(0),
 158 fAsym(pid.fAsym),
 159 fTPC(pid.fTPC),
 160 fTOF(pid.fTOF),
 161 fITS(pid.fITS),
 162 fTRD(pid.fTRD),
 163 fMatch(pid.fMatch),
 164 fForceTOFforKaons(pid.fForceTOFforKaons),
 165 fCompat(pid.fCompat),
 166 fPCompatTOF(pid.fPCompatTOF),
 167 fUseAsymTOF(pid.fUseAsymTOF),
 168 fLownSigmaTOF(pid.fLownSigmaTOF),
 169 fUpnSigmaTOF(pid.fUpnSigmaTOF),
 170 fLownSigmaCompatTOF(pid.fLownSigmaCompatTOF),
 171 fUpnSigmaCompatTOF(pid.fUpnSigmaCompatTOF),
 172 fnNSigmaCompat(pid.fnNSigmaCompat),
 173 fnSigmaCompat(0x0),
 174 fMC(pid.fMC),
 175 fOnePad(pid.fOnePad),
 176 fMCLowEn2011(pid.fMCLowEn2011),
 177 fppLowEn2011(pid.fppLowEn2011),
 178 fPbPb(pid.fPbPb),
 179 fTOFdecide(pid.fTOFdecide),
 180 fOldPid(pid.fOldPid),
 181 fPtThresholdTPC(pid.fPtThresholdTPC),
 182 fMaxTrackMomForCombinedPID(pid.fMaxTrackMomForCombinedPID),
 183 fPidResponse(0x0),
 184 fPidCombined(0x0),
 185 fTPCResponse(0x0),
 186 fCombDetectors(pid.fCombDetectors),
 187 fUseCombined(pid.fUseCombined),
 188 fDefaultPriors(pid.fDefaultPriors)
 189 {
 190
 191   fnSigmaCompat=new Double_t[fnNSigmaCompat];
 192   for(Int_t i=0;i<fnNSigmaCompat;i++){
 193     fnSigmaCompat[i]=pid.fnSigmaCompat[i];
 194   }
 195   fnSigma = new Double_t[fnNSigma];
 196   for(Int_t i=0;i<fnNSigma;i++){
 197     fnSigma[i]=pid.fnSigma[i];
 198   }
 199   fPriors = new Double_t[fnPriors];
 200   for(Int_t i=0;i<fnPriors;i++){
 201     fPriors[i]=pid.fPriors[i];
 202   }
 203   fPLimit = new Double_t[fnPLimit];
 204   for(Int_t i=0;i<fnPLimit;i++){
 205     fPLimit[i]=pid.fPLimit[i];
 206   }
 207   fPriors = new Double_t[fnPriors];
 208   for(Int_t i=0;i<fnPriors;i++){
 209     fPriors[i]=pid.fPriors[i];
 210   }
 211   for(Int_t i=0;i<AliPID::kSPECIES;i++){
 212     fPriorsH[i]=pid.fPriorsH[i];
 213   }
 214   for(Int_t i=0; i<3; i++){ // pi, K, proton
 215     fMaxnSigmaCombined[i]=pid.fMaxnSigmaCombined[i];
 216     fMinnSigmaTPC[i]=pid.fMinnSigmaTPC[i];
 217     fMaxnSigmaTPC[i]=pid.fMaxnSigmaTPC[i];
 218     fMinnSigmaTOF[i]=pid.fMinnSigmaTOF[i];
 219     fMaxnSigmaTOF[i]=pid.fMaxnSigmaTOF[i];
 220   }
 221
 222   //  if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
 223   fTPCResponse = new AliTPCPIDResponse();
 224   SetBetheBloch();
 225   fPidCombined = new AliPIDCombined();
 226   //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
 227   //fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
 228
 229   //Copy bands
 230   for (Int_t s=0;s<AliPID::kSPECIES;s++) {
 231     for (Int_t d=0;d<4;d++) {
 232       fIdBandMin[s][d]   = pid.fIdBandMin[s][d]   ? new TF1(*pid.fIdBandMin[s][d])   : NULL;
 233       fIdBandMax[s][d]   = pid.fIdBandMax[s][d]   ? new TF1(*pid.fIdBandMax[s][d])   : NULL;
 234       fCompBandMin[s][d] = pid.fCompBandMin[s][d] ? new TF1(*pid.fCompBandMin[s][d]) : NULL;
 235       fCompBandMax[s][d] = pid.fCompBandMax[s][d] ? new TF1(*pid.fCompBandMax[s][d]) : NULL;
 236     }
 237   }
 238
 239 }
 240 //----------------------
 241 Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
 242   // raw PID for single detectors, returns the particle type with smaller sigma
 243   Int_t specie=-1;
 244   if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
 245   if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
 246   if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
 247
 248   return specie;
 249
 250 }
 251 //---------------------------
 252 Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
 253   // checks if the track can be a kaon, raw PID applied for single detectors
 254   Int_t specie=0;
 255
 256   if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
 257   if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
 258   if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
 259
 260   if(specie==3) return kTRUE;
 261   return kFALSE;
 262 }
 263 //---------------------------
 264 Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
 265   // checks if the track can be a pion, raw PID applied for single detectors
 266
 267   Int_t specie=0;
 268
 269   if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
 270   if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
 271   if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
 272
 273   if(specie==2) return kTRUE;
 274   return kFALSE;
 275 }
 276 //---------------------------
 277 Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
 278   // checks if the track can be a proton raw PID applied for single detectors
 279
 280   Int_t specie=0;
 281   if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
 282   if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
 283   if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
 284
 285   if(specie==4) return kTRUE;
 286
 287   return kFALSE;
 288 }
 289 //--------------------------
 290 Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
 291   // checks if the track can be an electron raw PID applied for single detectors
 292
 293   Int_t specie=-1;
 294   if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
 295   if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
 296   if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
 297
 298   if(specie==0) return kTRUE;
 299
 300   return kFALSE;
 301 }
 302 //--------------------------
 303 Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
 304   // n-sigma cut, TPC PID
 305
 306   Double_t nsigma=-999.;
 307   Int_t pid=-1;
 308
 309   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
 310     Double_t nsigmaMin=999.;
 311     for(Int_t ipart=0;ipart<5;ipart++){
 312       if(GetnSigmaTPC(track,ipart,nsigma)==1){
 313         nsigma=TMath::Abs(nsigma);
 314         if((nsigma<nsigmaMin) && (nsigma<fnSigma[0])) {
 315           pid=ipart;
 316           nsigmaMin=nsigma;
 317         }
 318       }
 319     }
 320   }else{ // asks only for one particle specie
 321     if(GetnSigmaTPC(track,specie,nsigma)==1){
 322       nsigma=TMath::Abs(nsigma);
 323       if (nsigma>fnSigma[0]) pid=-1;
 324       else pid=specie;
 325     }
 326   }
 327
 328   return pid;
 329 }
 330 //----------------------------
 331 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
 332   // truncated mean, ITS PID
 333
 334   Double_t nsigma=-999.;
 335   Int_t pid=-1;
 336
 337   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
 338     Double_t nsigmaMin=999.;
 339     for(Int_t ipart=0;ipart<5;ipart++){
 340       if(GetnSigmaITS(track,ipart,nsigma)==1){
 341         nsigma=TMath::Abs(nsigma);
 342         if((nsigma<nsigmaMin) && (nsigma<fnSigma[4])) {
 343           pid=ipart;
 344           nsigmaMin=nsigma;
 345         }
 346       }
 347     }
 348   }else{ // asks only for one particle specie
 349     if(GetnSigmaITS(track,specie,nsigma)==1){
 350       nsigma=TMath::Abs(nsigma);
 351       if (nsigma>fnSigma[4]) pid=-1;
 352       else pid=specie;
 353     }
 354   }
 355
 356   return pid;
 357 }
 358 //----------------------------
 359 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
 360   // n-sigma cut, TOF PID
 361
 362   Double_t nsigma=-999.;
 363   Int_t pid=-1;
 364
 365   if(specie<0){
 366     Double_t nsigmaMin=999.;
 367     for(Int_t ipart=0;ipart<5;ipart++){
 368       if(GetnSigmaTOF(track,ipart,nsigma)==1){
 369         nsigma=TMath::Abs(nsigma);
 370         if((nsigma<nsigmaMin)&& (nsigma<fnSigma[3])){
 371           pid=ipart;
 372           nsigmaMin=nsigma;
 373         }
 374       }
 375     }
 376   }else{ // asks only for one particle specie
 377     Double_t nSigmaMin,nSigmaMax;
 378     if(fUseAsymTOF){
 379       nSigmaMin=fLownSigmaTOF;
 380       nSigmaMax=fUpnSigmaTOF;
 381     }else{
 382       nSigmaMin=-fnSigma[3];
 383       nSigmaMax=fnSigma[3];
 384     }
 385     if(GetnSigmaTOF(track,specie,nsigma)==1){
 386       if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
 387       else pid=specie;
 388     }
 389   }
 390   return pid;
 391 }
 392 //----------------------------
 393 Int_t AliAODPidHF::ApplyTOFCompatibilityBand(AliAODTrack *track,Int_t specie) const{
 394   // n-sigma cut, TOF PID
 395
 396   if(specie<0) return -1;
 397   Double_t nsigma=-999.;
 398   Int_t pid=-1;
 399
 400   Double_t nSigmaMin,nSigmaMax;
 401   if(fUseAsymTOF){
 402     nSigmaMin=fLownSigmaCompatTOF;
 403     nSigmaMax=fUpnSigmaCompatTOF;
 404   }else{
 405     nSigmaMin=-fnSigmaCompat[1];
 406     nSigmaMax=fnSigmaCompat[1];
 407   }
 408   if(GetnSigmaTOF(track,specie,nsigma)==1){
 409     if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
 410     else pid=specie;
 411   }
 412   return pid;
 413 }
 414 //------------------------------
 415 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
 416   // combined PID stored inside the AOD track
 417
 418   const Double_t *pid=track->PID();
 419   Float_t max=0.;
 420   Int_t k=-1;
 421   for (Int_t i=0; i<10; i++) {
 422     if (pid[i]>max) {k=i; max=pid[i];}
 423   }
 424
 425   if(k==2) type[0]=kTRUE;
 426   if(k==3) type[1]=kTRUE;
 427   if(k==4) type[2]=kTRUE;
 428
 429   return;
 430 }
 431 //--------------------------------
 432 Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
 433   // Check if the track is good for ITS PID
 434   AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kITS,track);
 435   if (status != AliPIDResponse::kDetPidOk) return kFALSE;
 436   return kTRUE;
 437 }
 438 //--------------------------------
 439 Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
 440   // Check if the track is good for TPC PID
 441   AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTPC,track);
 442   if (status != AliPIDResponse::kDetPidOk) return kFALSE;
 443   UInt_t nclsTPCPID = track->GetTPCsignalN();
 444   if(nclsTPCPID<fMinNClustersTPCPID) return kFALSE;
 445   return kTRUE;
 446 }
 447 //--------------------------------
 448 Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
 449   // Check if the track is good for TOF PID
 450   AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTOF,track);
 451   if (status != AliPIDResponse::kDetPidOk) return kFALSE;
 452   Float_t probMis = fPidResponse->GetTOFMismatchProbability(track);
 453   if (probMis > fCutTOFmismatch) return kFALSE;
 454   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0 &&
 455       track->GetStatus()&AliESDtrack::kITSrefit )   return kFALSE;
 456   return kTRUE;
 457 }
 458 //--------------------------------
 459 Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
 460   // Check if the track is good for TRD PID
 461   AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTRD,track);
 462   if (status != AliPIDResponse::kDetPidOk) return kFALSE;
 463   return kTRUE;
 464 }
 465 //--------------------------------
 466 Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
 467   // Quality cuts on the tracks, detector by detector
 468   if(detectors.Contains("ITS")) return CheckITSPIDStatus(track);
 469   else if(detectors.Contains("TPC")) return CheckTPCPIDStatus(track);
 470   else if(detectors.Contains("TOF")) return CheckTOFPIDStatus(track);
 471   else if(detectors.Contains("TRD")) return CheckTRDPIDStatus(track);
 472   else{
 473     AliError("Wrong detector name");
 474     return kFALSE;
 475   }
 476 }
 477 //--------------------------------------------
 478 Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
 479   // TPC nsigma cut PID, different sigmas in different p bins
 480
 481   AliAODPid *pidObj = track->GetDetPid();
 482   Double_t mom = pidObj->GetTPCmomentum();
 483   if(mom>fPtThresholdTPC) return kTRUE;
 484
 485   Double_t nsigma;
 486   if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
 487   nsigma=TMath::Abs(nsigma);
 488
 489
 490   if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
 491   if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
 492   if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
 493
 494   return kFALSE;
 495 }
 496 //------------------
 497 Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
 498   // combination of the PID info coming from TPC and TOF
 499
 500   Double_t ptrack=track->P();
 501   if(ptrack>fMaxTrackMomForCombinedPID) return 1;
 502
 503   Bool_t okTPC=CheckTPCPIDStatus(track);
 504   if(ptrack>fPtThresholdTPC) okTPC=kFALSE;
 505   Bool_t okTOF=CheckTOFPIDStatus(track);
 506
 507   if(fMatch==1){
 508     //TOF || TPC (a la' Andrea R.)
 509     // convention:
 510     // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
 511     // the method returns the sum of the response of the 2 detectors
 512
 513     if(fTPC && fTOF) {
 514       if(!okTPC && !okTOF) return 0;
 515     }
 516
 517     Int_t tTPCinfo=0;
 518     if(fTPC && okTPC){
 519       tTPCinfo=-1;
 520       if(fAsym) {
 521         if(TPCRawAsym(track,specie)) tTPCinfo=1;
 522       }else{
 523         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
 524       }
 525       if(fCompat && tTPCinfo<0){
 526         Double_t sig0tmp=fnSigma[0];
 527         SetSigma(0,fnSigmaCompat[0]);
 528         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
 529         SetSigma(0,sig0tmp);
 530       }
 531     }
 532
 533     Int_t tTOFinfo=0;
 534     if(fTOF){
 535       if(!okTOF && fTPC) return tTPCinfo;
 536       tTOFinfo=-1;
 537       if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
 538       if(fCompat && tTOFinfo>0){
 539         if(ptrack>fPCompatTOF) {
 540           if(ApplyTOFCompatibilityBand(track,specie)==specie) tTOFinfo=0;
 541         }
 542       }
 543     }
 544
 545
 546     if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
 547       if(!okTOF) return tTPCinfo;
 548       return tTOFinfo;
 549     }
 550
 551     if(tTPCinfo+tTOFinfo==0 && fITS){
 552       if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
 553       Int_t tITSinfo = -1;
 554       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 555       return tITSinfo;
 556     }
 557     return tTPCinfo+tTOFinfo;
 558   }
 559
 560   if(fMatch==2){
 561     //TPC & TOF (a la' Yifei)
 562     // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
 563     Int_t tTPCinfo=0;
 564
 565     if(fTPC && okTPC) {
 566       tTPCinfo=1;
 567       if(fAsym){
 568         if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
 569       }else{
 570         if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
 571       }
 572     }
 573
 574     Int_t tTOFinfo=1;
 575     if(fTOF){
 576       if(fTPC && !okTOF) return tTPCinfo;
 577       if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
 578     }
 579
 580     if(tTOFinfo==1 && tTPCinfo==1) return 1;
 581
 582     if(tTPCinfo+tTOFinfo==0 && fITS){
 583       if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
 584       Int_t tITSinfo = -1;
 585       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 586       return tITSinfo;
 587     }
 588     return -1;
 589   }
 590
 591   if(fMatch==3){
 592     //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
 593     // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
 594     if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
 595
 596
 597     Int_t tTPCinfo=-1;
 598     if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
 599       if(!okTPC) return 0;
 600       if(fAsym) {
 601         if(TPCRawAsym(track,specie)) tTPCinfo=1;
 602       }else{
 603         if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
 604       }
 605       return tTPCinfo;
 606     }
 607
 608     Int_t tTOFinfo=-1;
 609     if(ptrack>=fPLimit[1] && fTOF){
 610       if(!okTOF) return 0;
 611       if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
 612       return tTOFinfo;
 613     }
 614
 615     Int_t tITSinfo=-1;
 616     if(ptrack<fPLimit[0] && fITS){
 617       if(!CheckITSPIDStatus(track)) return 0;
 618       if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
 619       return tITSinfo;
 620     }
 621   }
 622
 623   if(fMatch==4 || fMatch==5){
 624
 625     // fMatch == 4 ---> "circular cut" in nSigmaTPC, nSimgaTOF plane
 626     //             ---> nsigmaTPC^2+nsigmaTOF^2 < cut^2
 627     // fMatch == 5 ---> "rectangular cut" in nSigmaTPC, nsigmaTOF plane
 628     //             ---> ns1<nSigmaTPC<NS1  && ns2<nSigmaTOF<NS2
 629
 630     Double_t nSigmaTPC=0.;
 631     if(okTPC) {
 632       nSigmaTPC=fPidResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)specie);
 633       if(nSigmaTPC<-990.) nSigmaTPC=0.;
 634     }
 635     Double_t nSigmaTOF=0.;
 636     if(okTOF) {
 637       nSigmaTOF=fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)specie);
 638     }
 639     Int_t iPart=specie-2; //species is 2 for pions,3 for kaons and 4 for protons
 640     if(iPart<0 || iPart>2) return -1;
 641     if(fMatch==4){
 642       Double_t nSigma2=nSigmaTPC*nSigmaTPC+nSigmaTOF*nSigmaTOF;
 643       if(nSigma2<fMaxnSigmaCombined[iPart]*fMaxnSigmaCombined[iPart]) return 1;
 644       else return -1;
 645     }
 646     else if(fMatch==5){
 647       if(fForceTOFforKaons && iPart==1 && !okTOF) return -1;
 648       if((nSigmaTPC>fMinnSigmaTPC[iPart] && nSigmaTPC<fMaxnSigmaTPC[iPart]) &&
 649          (nSigmaTOF>fMinnSigmaTOF[iPart] && nSigmaTOF<fMaxnSigmaTOF[iPart])) return 1;
 650       else return -1;
 651     }
 652   }
 653
 654   //Asymmetric cuts using user defined bands
 655   if (fMatch == 10) {
 656     if (fTPC && fTOF && !okTPC && !okTOF) {
 657       return 0;
 658     }
 659
 660     Int_t tTPCinfo = 0;
 661     if (fTPC && okTPC) {
 662       tTPCinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kTPC, track);
 663     }
 664
 665     Int_t tTOFinfo = 0;
 666     if (fTOF) {
 667       if (!okTOF && fTPC) {
 668         return tTPCinfo;
 669       }
 670       tTOFinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kTOF, track);
 671     }
 672
 673
 674     if (tTPCinfo+tTOFinfo == 0 && fTOFdecide) {
 675       if (!okTOF) {
 676         return tTPCinfo;
 677       }
 678       return tTOFinfo;
 679     }
 680
 681     if (tTPCinfo+tTOFinfo == 0 && fITS) {
 682       if (!CheckITSPIDStatus(track)) {
 683         return tTPCinfo+tTOFinfo;
 684       }
 685       Int_t tITSinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kITS, track);
 686       return tITSinfo;
 687     }
 688     return tTPCinfo+tTOFinfo;
 689   }
 690
 691   return -1;
 692
 693 }
 694 //----------------------------------
 695 Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
 696   // general method to compute PID
 697   if(fMatch>0){
 698     return MatchTPCTOF(track,specie);
 699   }else{
 700     if(fTPC && !fTOF && !fITS) {
 701       Int_t tTPCres=0;
 702       if(!fAsym){
 703         tTPCres=ApplyPidTPCRaw(track,specie);
 704         if(tTPCres==specie) return 1;
 705         else return tTPCres;
 706       }else{
 707         if(TPCRawAsym(track,specie)) tTPCres=1;
 708         else tTPCres=-1;
 709       }
 710       return tTPCres;
 711     }else if(fTOF && !fTPC && !fITS) {
 712       Int_t tTOFres=ApplyPidTOFRaw(track,specie);
 713       if(tTOFres==specie) return 1;
 714       else return tTOFres;
 715     }else if(fITS && !fTPC && !fTOF) {
 716       Int_t tITSres=ApplyPidITSRaw(track,specie);
 717       if(tITSres==specie) return 1;
 718       else return tITSres;
 719     }else{
 720       AliError("You should enable just one detector if you don't want to match");
 721       return 0;
 722     }
 723   }
 724 }
 725 //--------------------------------------------
 726 void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
 727   // TPC bethe bloch parameters
 728   if(fMC) {  // MC
 729
 730     if(fPbPb) { // PbPb MC
 731
 732       alephParameters[0] = 1.44405/50.;
 733       alephParameters[1] = 2.35409e+01;
 734       alephParameters[2] = TMath::Exp(-2.90330e+01);
 735       alephParameters[3] = 2.10681e+00;
 736       alephParameters[4] = 4.62254e+00;
 737
 738     } else {  // pp MC
 739       if(fMCLowEn2011){
 740         alephParameters[0]=0.0207667;
 741         alephParameters[1]=29.9936;
 742         alephParameters[2]=3.87866e-11;
 743         alephParameters[3]=2.17291;
 744         alephParameters[4]=7.1623;
 745       }else if(fOnePad){
 746         alephParameters[0]=0.029021;
 747         alephParameters[1]=25.4181;
 748         alephParameters[2]=4.66596e-08;
 749         alephParameters[3]=1.90008;
 750         alephParameters[4]=4.63783;
 751       }else{
 752         alephParameters[0] = 2.15898/50.;
 753         alephParameters[1] = 1.75295e+01;
 754         alephParameters[2] = 3.40030e-09;
 755         alephParameters[3] = 1.96178e+00;
 756         alephParameters[4] = 3.91720e+00;
 757       }
 758     }
 759
 760   } else { // Real Data
 761
 762     if(fOnePad) { // pp 1-pad (since LHC10d)
 763
 764       alephParameters[0] =1.34490e+00/50.;
 765       alephParameters[1] = 2.69455e+01;
 766       alephParameters[2] = TMath::Exp(-2.97552e+01);
 767       alephParameters[3] = 2.35339e+00;
 768       alephParameters[4] = 5.98079e+00;
 769
 770     } else if(fPbPb) { // PbPb
 771
 772       // alephParameters[0] = 1.25202/50.;
 773       // alephParameters[1] = 2.74992e+01;
 774       // alephParameters[2] = TMath::Exp(-3.31517e+01);
 775       // alephParameters[3] = 2.46246;
 776       // alephParameters[4] = 6.78938;
 777
 778       alephParameters[0] = 5.10207e+00/50.;
 779       alephParameters[1] = 7.94982e+00;
 780       alephParameters[2] = TMath::Exp(-9.07942e+00);
 781       alephParameters[3] = 2.38808e+00;
 782       alephParameters[4] = 1.68165e+00;
 783
 784     } else if(fppLowEn2011){ // pp low energy
 785
 786       alephParameters[0]=0.031642;
 787       alephParameters[1]=22.353;
 788       alephParameters[2]=4.16239e-12;
 789       alephParameters[3]=2.61952;
 790       alephParameters[4]=5.76086;
 791
 792     } else {  // pp no 1-pad (LHC10bc)
 793
 794       alephParameters[0] = 0.0283086/0.97;
 795       alephParameters[1] = 2.63394e+01;
 796       alephParameters[2] = 5.04114e-11;
 797       alephParameters[3] = 2.12543e+00;
 798       alephParameters[4] = 4.88663e+00;
 799
 800     }
 801
 802   }
 803
 804 }
 805
 806 //-----------------------
 807 void AliAODPidHF::SetBetheBloch() {
 808   // Set Bethe Bloch Parameters
 809
 810   Double_t alephParameters[5];
 811   GetTPCBetheBlochParams(alephParameters);
 812   fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 813
 814   return;
 815 }
 816
 817
 818 //--------------------------------------------------------------------------
 819 Int_t AliAODPidHF::GetnSigmaITS(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
 820   // get n sigma for ITS
 821
 822
 823   if (!CheckITSPIDStatus(track)) return -1;
 824
 825   Double_t nsigmaITS=-999;
 826
 827   if (fOldPid) {
 828     Double_t mom=track->P();
 829     AliAODPid *pidObj = track->GetDetPid();
 830     Double_t dedx=pidObj->GetITSsignal();
 831
 832     AliITSPIDResponse itsResponse;
 833     AliPID::EParticleType type=AliPID::EParticleType(species);
 834     nsigmaITS = itsResponse.GetNumberOfSigmas(mom,dedx,type);
 835
 836   } // old pid
 837   else { // new pid
 838
 839     AliPID::EParticleType type=AliPID::EParticleType(species);
 840     nsigmaITS = fPidResponse->NumberOfSigmasITS(track,type);
 841
 842   } //new pid
 843
 844   nsigma = nsigmaITS;
 845
 846   return 1;
 847
 848 }
 849 //--------------------------------------------------------------------------
 850 Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
 851   // get n sigma for TPC
 852
 853   if(!CheckTPCPIDStatus(track)) return -1;
 854
 855   Double_t nsigmaTPC=-999;
 856
 857   if(fOldPid){
 858     AliAODPid *pidObj = track->GetDetPid();
 859     Double_t dedx=pidObj->GetTPCsignal();
 860     Double_t mom = pidObj->GetTPCmomentum();
 861     if(mom>fPtThresholdTPC) return -2;
 862     UShort_t nTPCClus=pidObj->GetTPCsignalN();
 863     if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
 864     AliPID::EParticleType type=AliPID::EParticleType(species);
 865     nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
 866     nsigma=nsigmaTPC;
 867   } else{
 868     if(!fPidResponse) return -1;
 869     AliPID::EParticleType type=AliPID::EParticleType(species);
 870     nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
 871     nsigma=nsigmaTPC;
 872   }
 873   return 1;
 874 }
 875
 876 //-----------------------------
 877
 878 Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
 879   // get n sigma for TOF
 880
 881   if(!CheckTOFPIDStatus(track)) return -1;
 882
 883   if(fPidResponse){
 884     nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
 885     return 1;
 886   }else{
 887     AliFatal("To use TOF PID you need to attach AliPIDResponseTask");
 888     nsigma=-999.;
 889     return -1;
 890   }
 891 }
 892
 893 //-----------------------
 894 Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
 895   // Exclude a given hypothesis (labelTracks) in detector
 896
 897   if (detectors.Contains("ITS")) {
 898
 899     AliInfo("Nothing to be done");
 900     /*
 901      Double_t nsigma=0.;
 902      if (GetnSigmaITS(track,labelTrack,nsigma)==1){
 903      if(nsigma>nsigmaCut) return kTRUE;
 904      }
 905      */
 906     return kFALSE;
 907
 908   } else if (detectors.Contains("TPC")) {
 909
 910     Double_t nsigma=0.;
 911     if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
 912       if(nsigma>nsigmaCut) return kTRUE;
 913     }
 914     return kFALSE;
 915
 916   } else if (detectors.Contains("TOF")) {
 917
 918     Double_t nsigma=0.;
 919     if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
 920       if(nsigma>nsigmaCut) return kTRUE;
 921     }
 922     return kFALSE;
 923
 924   }
 925   return kFALSE;
 926
 927 }
 928 //-----------------------
 929 Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
 930   // TOF proton compatibility
 931
 932   if(!CheckTOFPIDStatus(track)) return 0;
 933
 934   Double_t nsigma;
 935   if(GetnSigmaTOF(track,3,nsigma)==1){
 936     if(nsigma>nsigmaK) return kTRUE;
 937   }
 938   return kFALSE;
 939   /*  Double_t time[AliPID::kSPECIESN];
 940    Double_t sigmaTOFPid[AliPID::kSPECIES];
 941    AliAODPid *pidObj = track->GetDetPid();
 942    pidObj->GetIntegratedTimes(time);
 943    Double_t sigTOF=pidObj->GetTOFsignal();
 944
 945    AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
 946    if (event) {
 947    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
 948    if (tofH && fPidResponse) {
 949    AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
 950    sigTOF -= TOFres.GetStartTime(track->P());
 951    sigmaTOFPid[3]=TOFres.GetExpectedSigma(track->P(),time[3],AliPID::ParticleMass(3));
 952    }
 953    else  pidObj->GetTOFpidResolution(sigmaTOFPid);
 954    } else  pidObj->GetTOFpidResolution(sigmaTOFPid);
 955    Double_t sigmaTOFtrack;
 956    if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
 957    else sigmaTOFtrack=fTOFSigma;  // backward compatibility for old AODs
 958
 959    if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
 960
 961    return kFALSE;
 962    */
 963 }
 964
 965 //--------------------------------------------------------------------------
 966 void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
 967
 968         //
 969         // method setting the prior distributions to the AliPIDCombined object of the AliAODPidHF data member
 970         // all the checks are done directly in the AliPIDCombined object
 971         //
 972
 973         GetPidCombined()->SetPriorDistribution(type,prior);
 974 }
 975 //--------------------------------------------------------------------------
 976 void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
 977
 978         //
 979         // Drawing prior distribution for type "type"
 980
 981         new TCanvas();
 982         GetPidCombined()->GetPriorDistribution(type)->Draw();
 983 }
 984
 985 //-----------------------------
 986 void AliAODPidHF::SetPriorsHistos(TString priorFileName){
 987   // Set histograms with priors
 988
 989   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
 990     if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
 991     TString nt ="name";
 992     nt+="_prior_";
 993     nt+=AliPID::ParticleName(ispecies);
 994   }
 995   TDirectory *current = gDirectory;
 996   TFile *priorFile=TFile::Open(priorFileName);
 997   if (priorFile) {
 998     TH1F* h3=static_cast<TH1F*>(priorFile->Get("priors3step9"));
 999     TH1F* h2=static_cast<TH1F*>(priorFile->Get("priors2step9"));
1000     TH1F* h1=static_cast<TH1F*>(priorFile->Get("priors1step9"));
1001     current->cd();
1002     fPriorsH[AliPID::kProton] = new TH1F(*h3);
1003     fPriorsH[AliPID::kKaon  ] = new TH1F(*h2);
1004     fPriorsH[AliPID::kPion  ] = new TH1F(*h1);
1005     priorFile->Close();
1006     delete priorFile;
1007     TF1 *salt=new TF1("salt","1.e-10",0,10);
1008     fPriorsH[AliPID::kProton]->Add(salt);
1009     fPriorsH[AliPID::kKaon  ]->Add(salt);
1010     fPriorsH[AliPID::kPion  ]->Add(salt);
1011     delete salt;
1012   }
1013 }
1014 //----------------------------------
1015 void AliAODPidHF::SetUpCombinedPID(){
1016   // Configuration of combined Bayesian PID
1017
1018   fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
1019   if(!fDefaultPriors){
1020         for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
1021         fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
1022         }
1023   }else{
1024         fPidCombined->SetDefaultTPCPriors();
1025   }
1026   switch (fCombDetectors){
1027     case kTPCTOF:
1028       fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1029       break;
1030     case kTPCITS:
1031       fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
1032       break;
1033     case kTPC:
1034       fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1035       break;
1036     case kTOF:
1037       fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1038       break;
1039   }
1040 }
1041
1042
1043 //-----------------------------
1044 void AliAODPidHF::PrintAll() const {
1045   // print the configuration
1046   printf("Detectors used for PID: ");
1047   if(fITS) printf("ITS ");
1048   if(fTPC) printf("TPC ");
1049   if(fTRD) printf("TRD ");
1050   if(fTOF) printf("TOF ");
1051   printf("\n");
1052   printf("Minimum TPC PID clusters = %d\n",fMinNClustersTPCPID);
1053   printf("Maximum momentum for using TPC PID = %f\n",fPtThresholdTPC);
1054   printf("TOF Mismatch probablility cut = %f\n",fCutTOFmismatch);
1055   printf("Maximum momentum for combined PID TPC PID = %f\n",fMaxTrackMomForCombinedPID);
1056   if(fOldPid){
1057     printf("Use OLD PID");
1058     printf("  fMC = %d\n",fMC);
1059     printf("  fPbPb = %d\n",fPbPb);
1060     printf("  fOnePad = %d\n",fOnePad);
1061     printf("  fMCLowEn2011 = %d\n",fMCLowEn2011);
1062     printf("  fppLowEn2011 = %d\n",fppLowEn2011);
1063   }
1064   printf("--- Matching algorithm = %d ---\n",fMatch);
1065   if(fMatch==1){
1066     if(fITS) printf("nSigmaITS = %.2f\n",fnSigma[4]);
1067     if(fTOF){
1068       printf("nSigmaTOF = %.2f\n",fnSigma[3]);
1069       if(fCompat) printf("Compatibility band at nSigmaTOF=%.2f for p>%.2f\n",fnSigmaCompat[1],fPCompatTOF);
1070     }
1071     if(fTPC){
1072       if(fAsym){
1073         printf("nSigmaTPC:\n");
1074         printf("   pt<%.2f      \t nsigmaTPC= %.2f\n",fPLimit[0],fnSigma[0]);
1075         printf("   %.2f<pt<%.2f \t nsigmaTPC= %.2f\n",fPLimit[0],fPLimit[1],fnSigma[1]);
1076         printf("   pt>%.2f      \t nsigmaTPC= %.2f\n",fPLimit[1],fnSigma[2]);
1077       }else{
1078         printf("nSigmaTPC = %.2f\n",fnSigma[0]);
1079       }
1080       if(fCompat) printf("Compatibility band at nSigmaTPC=%.2f\n",fnSigmaCompat[0]);
1081     }
1082   }else if(fMatch==4){
1083     printf("Cuts on sqrt(nSigmaTPC^2+nSigmaTOF^2):\n");
1084     printf(" Pions:   nSigma = %.2f\n",fMaxnSigmaCombined[0]);
1085     printf(" Kaons:   nSigma = %.2f\n",fMaxnSigmaCombined[1]);
1086     printf(" Protons: nSigma = %.2f\n",fMaxnSigmaCombined[2]);
1087   }else if(fMatch==5){
1088     printf("nSigma ranges:\n");
1089     printf(" Pions:   %.2f<nSigmaTPC<%.2f   %.2f<nSigmaTOF<%.2f\n",
1090            fMinnSigmaTPC[0],fMaxnSigmaTPC[0],fMinnSigmaTOF[0],fMaxnSigmaTOF[0]);
1091     printf(" Kaons:   %.2f<nSigmaTPC<%.2f   %.2f<nSigmaTOF<%.2f\n",
1092            fMinnSigmaTPC[1],fMaxnSigmaTPC[1],fMinnSigmaTOF[1],fMaxnSigmaTOF[1]);
1093     printf(" Protons: %.2f<nSigmaTPC<%.2f   %.2f<nSigmaTOF<%.2f\n",
1094            fMinnSigmaTPC[2],fMaxnSigmaTPC[2],fMinnSigmaTOF[2],fMaxnSigmaTOF[2]);
1095   } else if (fMatch == 10) {
1096     printf("Asymmetric PID using identification/compatibility bands as a function of track momentum p\n");
1097     printf("The following bands are set:\n");
1098     TString species[] = {"electron", "muon", "pion", "kaon", "proton"};
1099     TString detectors[] = {"ITS", "TPC", "TRD", "TOF"};
1100     for (Int_t s=0;s<AliPID::kSPECIES;s++) {
1101       for (Int_t d=0;d<4;d++) {
1102         if (fIdBandMin[s][d] && fIdBandMax[s][d]) {
1103           printf("  Identification band %s %s\n", species[s].Data(), detectors[d].Data());
1104         }
1105         if (fCompBandMin[s][d] && fCompBandMax[s][d]) {
1106           printf("  Compatibility band %s %s\n", species[s].Data(), detectors[d].Data());
1107         }
1108       }
1109     }
1110   }
1111 }
1112
1113 //------------------
1114 void AliAODPidHF::SetIdBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max) {
1115   int spe = (int) specie;
1116   int det = (int) detector;
1117
1118   if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1119     AliError("Identification band not set");
1120     return;
1121   }
1122
1123   TAxis *axis;
1124   HistFunc *histFunc;
1125
1126   axis = min->GetXaxis();
1127   histFunc = new HistFunc(min);
1128   TF1 *minFunc = new TF1(Form("IdMin_%d_%d", spe, det), histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1129
1130   axis = max->GetXaxis();
1131   histFunc = new HistFunc(max);
1132   TF1 *maxFunc = new TF1(Form("IdMax_%d_%d", spe, det), histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1133
1134   SetIdBand(specie, detector, minFunc, maxFunc);
1135 }
1136
1137 //------------------
1138 void AliAODPidHF::SetIdBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TF1 *min, TF1 *max) {
1139   int spe = (int) specie;
1140   int det = (int) detector;
1141
1142   if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1143     AliError("Identification band not set");
1144     return;
1145   }
1146
1147   if (fIdBandMin[spe][det]) {
1148     delete fIdBandMin[spe][det];
1149   }
1150   fIdBandMin[spe][det] = new TF1(*min);
1151
1152   if (fIdBandMax[spe][det]) {
1153     delete fIdBandMax[spe][det];
1154   }
1155   fIdBandMax[spe][det] = new TF1(*max);
1156 }
1157
1158 //------------------
1159 void AliAODPidHF::SetCompBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max) {
1160   int spe = (int) specie;
1161   int det = (int) detector;
1162
1163   if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1164     AliError("Compatibility band not set");
1165     return;
1166   }
1167
1168   TAxis *axis;
1169   HistFunc *histFunc;
1170
1171   axis = min->GetXaxis();
1172   histFunc = new HistFunc(min);
1173   TF1 *minFunc = new TF1(Form("CompMin_%d_%d", spe, det), histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1174
1175   axis = max->GetXaxis();
1176   histFunc = new HistFunc(max);
1177   TF1 *maxFunc = new TF1(Form("CompMax_%d_%d", spe, det), histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1178
1179   SetCompBand(specie, detector, minFunc, maxFunc);
1180 }
1181
1182 //------------------
1183 void AliAODPidHF::SetCompBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TF1 *min, TF1 *max) {
1184   int spe = (int) specie;
1185   int det = (int) detector;
1186
1187   if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1188     AliError("Compatibility band not set");
1189     return;
1190   }
1191
1192   if (fCompBandMin[spe][det]) {
1193     delete fCompBandMin[spe][det];
1194   }
1195   fCompBandMin[spe][det] = new TF1(*min);
1196
1197   if (fCompBandMax[spe][det]) {
1198     delete fCompBandMax[spe][det];
1199   }
1200   fCompBandMax[spe][det] = new TF1(*max);
1201 }
1202
1203 //------------------
1204 Bool_t AliAODPidHF::CheckDetectorPIDStatus(AliPIDResponse::EDetector detector, AliAODTrack* track) {
1205   switch (detector) {
1206     case AliPIDResponse::kITS:
1207       return CheckITSPIDStatus(track);
1208       break;
1209     case AliPIDResponse::kTPC:
1210       return CheckTPCPIDStatus(track);
1211       break;
1212     case AliPIDResponse::kTRD:
1213       return CheckTRDPIDStatus(track);
1214       break;
1215     case AliPIDResponse::kTOF:
1216       return CheckTOFPIDStatus(track);
1217       break;
1218     default:
1219       return kFALSE;
1220       break;
1221   }
1222 }
1223
1224 //------------------
1225 Float_t AliAODPidHF::NumberOfSigmas(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track) {
1226   switch (detector) {
1227     case AliPIDResponse::kITS:
1228       return fPidResponse->NumberOfSigmasITS(track, specie);
1229       break;
1230     case AliPIDResponse::kTPC:
1231       return fPidResponse->NumberOfSigmasTPC(track, specie);
1232       break;
1233     case AliPIDResponse::kTOF:
1234       return fPidResponse->NumberOfSigmasTOF(track, specie);
1235       break;
1236     default:
1237       return -999.;
1238       break;
1239   }
1240 }
1241
1242 //------------------
1243 Int_t AliAODPidHF::CheckBands(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track) {
1244   //Return: -1 for no match, 0 for comp. band and 1 for id band
1245
1246   Int_t spe = (Int_t) specie;
1247   Int_t det = (Int_t) detector;
1248
1249   if (!fPidResponse || spe >= AliPID::kSPECIES) {
1250     return -1;
1251   }
1252
1253   if (!CheckDetectorPIDStatus(detector, track)) {
1254     return 1;
1255   }
1256
1257   double P = track->P();
1258
1259   Float_t nSigma = NumberOfSigmas(specie, detector, track);
1260   Float_t minContent, maxContent;
1261
1262   //Check if within identification band, return 1
1263   TF1 *IdBandMin = fIdBandMin[spe][det];
1264   TF1 *IdBandMax = fIdBandMax[spe][det];
1265
1266   if (IdBandMin && IdBandMax) {
1267     minContent = IdBandMin->Eval(P);
1268     maxContent = IdBandMax->Eval(P);
1269     if ((minContent == 0 || nSigma >= minContent) && (maxContent == 0 || nSigma <= maxContent)) {
1270       return 1;
1271     }
1272   }
1273
1274   //Check if within compatibility band, return 0
1275   TF1 *CompBandMin = fCompBandMin[spe][det];
1276   TF1 *CompBandMax = fCompBandMax[spe][det];
1277
1278   if (CompBandMin && CompBandMax) {
1279     minContent = CompBandMin->Eval(P);
1280     maxContent = CompBandMax->Eval(P);
1281     if ((minContent == 0 || nSigma >= minContent) && (maxContent == 0 || nSigma <= maxContent)) {
1282       return 0;
1283     }
1284   }
1285
1286   //No bands set: don't perform PID
1287   if (!IdBandMin && !IdBandMax && !CompBandMin && !CompBandMax) {
1288     AliWarning(Form("No identification & compatibility band set for specie=%d detector=%d", spe, det));
1289     return 1;
1290   }
1291
1292   //Otherwise, return -1
1293   return -1;
1294 }
1295
1296 //------------------
1297 void AliAODPidHF::SetShiftedAsymmetricPID() {
1298   SetMatch(10);
1299   SetTPC(kTRUE);
1300   SetTOF(kTRUE);
1301
1302   //TPC K: shift by -0.2
1303   TF1 *TPCCompBandMinK = new TF1("TPCCompBandMinK", "[0]", 0.3, 4); TPCCompBandMinK->SetParameter(0, -3.2);
1304   TF1 *TPCCompBandMaxK = new TF1("TPCCompBandMaxK", "[0]", 0.3, 4); TPCCompBandMaxK->SetParameter(0, 2.8);
1305   SetCompBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCCompBandMinK, TPCCompBandMaxK);
1306
1307   TF1 *TPCIdBandMinK = new TF1("TPCIdBandMinK", "[0]", 0.3, 4); TPCIdBandMinK->SetParameter(0, -2.2);
1308   TF1 *TPCIdBandMaxK = new TF1("TPCIdBandMaxK", "[0]", 0.3, 4); TPCIdBandMaxK->SetParameter(0, 1.8);
1309   SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1310
1311   //TPC pi: shift by -0.14
1312   TF1 *TPCCompBandMinPi = new TF1("TPCCompBandMinPi", "[0]", 0.3, 4); TPCCompBandMinPi->SetParameter(0, -3.14);
1313   TF1 *TPCCompBandMaxPi = new TF1("TPCCompBandMaxPi", "[0]", 0.3, 4); TPCCompBandMaxPi->SetParameter(0, 2.86);
1314   SetCompBand(AliPID::kPion, AliPIDResponse::kTPC, TPCCompBandMinPi, TPCCompBandMaxPi);
1315
1316   TF1 *TPCIdBandMinPi = new TF1("TPCIdBandMinPi", "[0]", 0.3, 4); TPCIdBandMinPi->SetParameter(0, -2.14);
1317   TF1 *TPCIdBandMaxPi = new TF1("TPCIdBandMaxPi", "[0]", 0.3, 4); TPCIdBandMaxPi->SetParameter(0, 1.86);
1318   SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinPi, TPCIdBandMaxPi);
1319
1320   //TOF K: shift by -0.1
1321   TF1 *TOFCompBandMinK = new TF1("TOFCompBandMinK", "[0]", 2, 50); TOFCompBandMinK->SetParameter(0, -3.1);
1322   TF1 *TOFCompBandMaxK = new TF1("TOFCompBandMaxK", "[0]", 2, 50); TOFCompBandMaxK->SetParameter(0, 2.9);
1323   SetCompBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFCompBandMinK, TOFCompBandMaxK);
1324
1325   TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0.3, 2); TOFIdBandMinK->SetParameter(0, -3.1);
1326   TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0.3, 2); TOFIdBandMaxK->SetParameter(0, 2.9);
1327   SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1328
1329   //TOF pi: shift by -0.15
1330   TF1 *TOFCompBandMinPi = new TF1("TOFCompBandMinPi", "[0]", 2, 50); TOFCompBandMinPi->SetParameter(0, -3.15);
1331   TF1 *TOFCompBandMaxPi = new TF1("TOFCompBandMaxPi", "[0]", 2, 50); TOFCompBandMaxPi->SetParameter(0, 2.85);
1332   SetCompBand(AliPID::kPion, AliPIDResponse::kTOF, TOFCompBandMinPi, TOFCompBandMaxPi);
1333
1334   TF1 *TOFIdBandMinPi = new TF1("TOFIdBandMinPi", "[0]", 0.3, 2); TOFIdBandMinPi->SetParameter(0, -3.15);
1335   TF1 *TOFIdBandMaxPi = new TF1("TOFIdBandMaxPi", "[0]", 0.3, 2); TOFIdBandMaxPi->SetParameter(0, 2.85);
1336   SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinPi, TOFIdBandMaxPi);
1337 }
1338
1339 //------------------
1340 void AliAODPidHF::SetIdAsymmetricPID() {
1341   //Set identification bands
1342
1343   SetMatch(10);
1344   SetTPC(kTRUE);
1345   SetTOF(kTRUE);
1346
1347   //TPC K
1348   Double_t TPCIdBandMinKBins[] = {0.3, 0.4, 0.5, 0.6, 0.9, 4};
1349   TH1F *TPCIdBandMinK = new TH1F("TPCIdBandMinK", "TPC Id Band Min K", 5, TPCIdBandMinKBins);
1350   TPCIdBandMinK->SetBinContent(1, -3); //0.3-0.4
1351   TPCIdBandMinK->SetBinContent(2, -2); //0.4-0.5
1352   TPCIdBandMinK->SetBinContent(3, -3); //0.5-0.6
1353   TPCIdBandMinK->SetBinContent(4, -2); //0.6-0.9
1354   TPCIdBandMinK->SetBinContent(5, -3); //0.9-4
1355
1356   Double_t TPCIdBandMaxKBins[] = {0.3, 0.6, 0.7, 4};
1357   TH1F *TPCIdBandMaxK = new TH1F("TPCIdBandMaxK", "TPC Id Band Max K", 3, TPCIdBandMaxKBins);
1358   TPCIdBandMaxK->SetBinContent(1, 3); //0.3-0.6
1359   TPCIdBandMaxK->SetBinContent(2, 2); //0.6-0.7
1360   TPCIdBandMaxK->SetBinContent(3, 3); //0.7-4
1361
1362   SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1363   GetIdBandMin(AliPID::kKaon, AliPIDResponse::kTPC)->SetNpx(1000);
1364   GetIdBandMax(AliPID::kKaon, AliPIDResponse::kTPC)->SetNpx(1000);
1365
1366   //TPC pi
1367   Double_t TPCIdBandMinpiBins[] = {0.3, 4};
1368   TH1F *TPCIdBandMinpi = new TH1F("TPCIdBandMinpi", "TPC Id Band Min pi", 1, TPCIdBandMinpiBins);
1369   TPCIdBandMinpi->SetBinContent(1, -3); //0.3-4
1370
1371   Double_t TPCIdBandMaxpiBins[] = {0.3, 0.7, 0.9, 1.3, 1.4, 4};
1372   TH1F *TPCIdBandMaxpi = new TH1F("TPCIdBandMaxpi", "TPC Id Band Max pi", 5, TPCIdBandMaxpiBins);
1373   TPCIdBandMaxpi->SetBinContent(1, 3); //0.3-0.7
1374   TPCIdBandMaxpi->SetBinContent(2, 2); //0.7-0.9
1375   TPCIdBandMaxpi->SetBinContent(3, 3); //0.9-1.3
1376   TPCIdBandMaxpi->SetBinContent(4, 2); //1.3-1.4
1377   TPCIdBandMaxpi->SetBinContent(5, 3); //1.4-4
1378
1379   SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinpi, TPCIdBandMaxpi);
1380   GetIdBandMin(AliPID::kPion, AliPIDResponse::kTPC)->SetNpx(1000);
1381   GetIdBandMax(AliPID::kPion, AliPIDResponse::kTPC)->SetNpx(1000);
1382
1383   //TOF K
1384   TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0, 50); TOFIdBandMinK->SetParameter(0, -3);
1385   TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0, 50); TOFIdBandMaxK->SetParameter(0, 3);
1386
1387   SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1388
1389   //TOF pi
1390   TF1 *TOFIdBandMinPi = new TF1("TOFIdBandMinPi", "[0]", 0, 50); TOFIdBandMinPi->SetParameter(0, -3);
1391   TF1 *TOFIdBandMaxPi = new TF1("TOFIdBandMaxPi", "[0]", 0, 50); TOFIdBandMaxPi->SetParameter(0, 3);
1392
1393   SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinPi, TOFIdBandMaxPi);
1394 }
1395
1396 //------------------
1397 void AliAODPidHF::SetIdCompAsymmetricPID() {
1398   //Set compatibility and identification bands
1399
1400   SetMatch(10);
1401   SetTPC(kTRUE);
1402   SetTOF(kTRUE);
1403
1404   //TPC K
1405   TF1 *TPCCompBandMinK = new TF1("TPCCompBandMinK", "[0]", 0.3, 4); TPCCompBandMinK->SetParameter(0, -3);
1406   TF1 *TPCCompBandMaxK = new TF1("TPCCompBandMaxK", "[0]", 0.3, 4); TPCCompBandMaxK->SetParameter(0, 3);
1407
1408   SetCompBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCCompBandMinK, TPCCompBandMaxK);
1409
1410   Double_t TPCIdBandMinKBins[6] = {0.3, 0.45, 0.55, 0.7, 1.1, 4};
1411   TH1F *TPCIdBandMinK = new TH1F("TPCIdBandMinK", "TPC Id Band Min K", 5, TPCIdBandMinKBins);
1412   TPCIdBandMinK->SetBinContent(1, -2); //0.3 -0.45
1413   TPCIdBandMinK->SetBinContent(2, -1); //0.45-0.55
1414   TPCIdBandMinK->SetBinContent(3, -2); //0.55-0.7
1415   TPCIdBandMinK->SetBinContent(4, -1); //0.7 -1.1
1416   TPCIdBandMinK->SetBinContent(5, -2); //1.1 -4
1417
1418   Double_t TPCIdBandMaxKBins[4] = {0.3, 0.5, 0.7, 4};
1419   TH1F *TPCIdBandMaxK = new TH1F("TPCIdBandMaxK", "TPC Id Band Max K", 3, TPCIdBandMaxKBins);
1420   TPCIdBandMaxK->SetBinContent(1, 2); //0.3-0.5
1421   TPCIdBandMaxK->SetBinContent(2, 1); //0.5-0.7
1422   TPCIdBandMaxK->SetBinContent(3, 2); //0.7-4
1423
1424   SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1425   GetIdBandMin(AliPID::kKaon, AliPIDResponse::kTPC)->SetNpx(1000);
1426   GetIdBandMax(AliPID::kKaon, AliPIDResponse::kTPC)->SetNpx(1000);
1427
1428   //TPC pi
1429   TF1 *TPCCompBandMinpi = new TF1("TPCCompBandMinpi", "[0]", 0.3, 4); TPCCompBandMinpi->SetParameter(0, -3);
1430   TF1 *TPCCompBandMaxpi = new TF1("TPCCompBandMaxpi", "[0]", 0.3, 4); TPCCompBandMaxpi->SetParameter(0, 3);
1431
1432   SetCompBand(AliPID::kPion, AliPIDResponse::kTPC, TPCCompBandMinpi, TPCCompBandMaxpi);
1433
1434   Double_t TPCIdBandMinpiBins[2] = {0.3, 4};
1435   TH1F *TPCIdBandMinpi = new TH1F("TPCIdBandMinpi", "TPC Id Band Min pi", 1, TPCIdBandMinpiBins);
1436   TPCIdBandMinpi->SetBinContent(1, -2); //0.3-4
1437
1438   Double_t TPCIdBandMaxpiBins[4] = {0.3, 0.7, 1.7, 4};
1439   TH1F *TPCIdBandMaxpi = new TH1F("TPCIdBandMaxpi", "TPC Id Band Max pi", 3, TPCIdBandMaxpiBins);
1440   TPCIdBandMaxpi->SetBinContent(1, 2); //0.3-0.7
1441   TPCIdBandMaxpi->SetBinContent(2, 1); //0.7-1.7
1442   TPCIdBandMaxpi->SetBinContent(3, 2); //1.7-4
1443
1444   SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinpi, TPCIdBandMaxpi);
1445   GetIdBandMin(AliPID::kPion, AliPIDResponse::kTPC)->SetNpx(1000);
1446   GetIdBandMax(AliPID::kPion, AliPIDResponse::kTPC)->SetNpx(1000);
1447
1448   //TOF K
1449   TF1 *TOFCompBandMinK = new TF1("TOFCompBandMinK", "[0]", 2, 50); TOFCompBandMinK->SetParameter(0, -3);
1450   TF1 *TOFCompBandMaxK = new TF1("TOFCompBandMaxK", "[0]", 2, 50); TOFCompBandMaxK->SetParameter(0, 3);
1451
1452   SetCompBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFCompBandMinK, TOFCompBandMaxK);
1453
1454   TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0.3, 2); TOFIdBandMinK->SetParameter(0, -3);
1455   TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0.3, 2); TOFIdBandMaxK->SetParameter(0, 3);
1456
1457   SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1458
1459   //TOF pi
1460   TF1 *TOFCompBandMinpi = new TF1("TOFCompBandMinpi", "[0]", 2, 50); TOFCompBandMinpi->SetParameter(0, -3);
1461   TF1 *TOFCompBandMaxpi = new TF1("TOFCompBandMaxpi", "[0]", 2, 50); TOFCompBandMaxpi->SetParameter(0, 3);
1462
1463   SetCompBand(AliPID::kPion, AliPIDResponse::kTOF, TOFCompBandMinpi, TOFCompBandMaxpi);
1464
1465   TF1 *TOFIdBandMinpi = new TF1("TOFIdBandMinpi", "[0]", 0.3, 2); TOFIdBandMinpi->SetParameter(0, -3);
1466   TF1 *TOFIdBandMaxpi = new TF1("TOFIdBandMaxpi", "[0]", 0.3, 2); TOFIdBandMaxpi->SetParameter(0, 3);
1467
1468   SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinpi, TOFIdBandMaxpi);
1469 }