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