STEER/STEERBase/AliPIDResponse.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
  17
  18 //-----------------------------------------------------------------
  19 //        Base class for handling the pid response               //
  20 //        functions of all detectors                             //
  21 //        and give access to the nsigmas                         //
  22 //                                                               //
  23 //   Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch //
  24 //-----------------------------------------------------------------
  25
  26 #include <TList.h>
  27 #include <TObjArray.h>
  28 #include <TPRegexp.h>
  29 #include <TF1.h>
  30 #include <TSpline.h>
  31 #include <TFile.h>
  32 #include <TArrayI.h>
  33 #include <TArrayF.h>
  34
  35 #include <AliVEvent.h>
  36 #include <AliVTrack.h>
  37 #include <AliLog.h>
  38 #include <AliPID.h>
  39 #include <AliOADBContainer.h>
  40 #include <AliTRDPIDResponseObject.h>
  41 #include <AliTOFPIDParams.h>
  42
  43 #include "AliPIDResponse.h"
  44 #include "AliDetectorPID.h"
  45
  46 #include "AliCentrality.h"
  47
  48 ClassImp(AliPIDResponse);
  49
  50 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
  51 TNamed("PIDResponse","PIDResponse"),
  52 fITSResponse(isMC),
  53 fTPCResponse(),
  54 fTRDResponse(),
  55 fTOFResponse(),
  56 fEMCALResponse(),
  57 fRange(5.),
  58 fITSPIDmethod(kITSTruncMean),
  59 fIsMC(isMC),
  60 fOADBPath(),
  61 fCustomTPCpidResponse(),
  62 fBeamType("PP"),
  63 fLHCperiod(),
  64 fMCperiodTPC(),
  65 fMCperiodUser(),
  66 fCurrentFile(),
  67 fRecoPass(0),
  68 fRecoPassUser(-1),
  69 fRun(0),
  70 fOldRun(0),
  71 fArrPidResponseMaster(0x0),
  72 fResolutionCorrection(0x0),
  73 fTRDPIDResponseObject(0x0),
  74 fTOFtail(1.1),
  75 fTOFPIDParams(0x0),
  76 fEMCALPIDParams(0x0),
  77 fCurrentEvent(0x0),
  78 fCurrCentrality(0.0),
  79 fTuneMConData(kFALSE)
  80 {
  81   //
  82   // default ctor
  83   //
  84   AliLog::SetClassDebugLevel("AliPIDResponse",0);
  85   AliLog::SetClassDebugLevel("AliESDpid",0);
  86   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
  87
  88   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
  89 }
  90
  91 //______________________________________________________________________________
  92 AliPIDResponse::~AliPIDResponse()
  93 {
  94   //
  95   // dtor
  96   //
  97   delete fArrPidResponseMaster;
  98   delete fTRDPIDResponseObject;
  99   delete fTOFPIDParams;
 100 }
 101
 102 //______________________________________________________________________________
 103 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
 104 TNamed(other),
 105 fITSResponse(other.fITSResponse),
 106 fTPCResponse(other.fTPCResponse),
 107 fTRDResponse(other.fTRDResponse),
 108 fTOFResponse(other.fTOFResponse),
 109 fEMCALResponse(other.fEMCALResponse),
 110 fRange(other.fRange),
 111 fITSPIDmethod(other.fITSPIDmethod),
 112 fIsMC(other.fIsMC),
 113 fOADBPath(other.fOADBPath),
 114 fCustomTPCpidResponse(other.fCustomTPCpidResponse),
 115 fBeamType("PP"),
 116 fLHCperiod(),
 117 fMCperiodTPC(),
 118 fMCperiodUser(other.fMCperiodUser),
 119 fCurrentFile(),
 120 fRecoPass(0),
 121 fRecoPassUser(other.fRecoPassUser),
 122 fRun(0),
 123 fOldRun(0),
 124 fArrPidResponseMaster(0x0),
 125 fResolutionCorrection(0x0),
 126 fTRDPIDResponseObject(0x0),
 127 fTOFtail(1.1),
 128 fTOFPIDParams(0x0),
 129 fEMCALPIDParams(0x0),
 130 fCurrentEvent(0x0),
 131 fCurrCentrality(0.0),
 132 fTuneMConData(kFALSE)
 133 {
 134   //
 135   // copy ctor
 136   //
 137   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 138 }
 139
 140 //______________________________________________________________________________
 141 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
 142 {
 143   //
 144   // copy ctor
 145   //
 146   if(this!=&other) {
 147     delete fArrPidResponseMaster;
 148     TNamed::operator=(other);
 149     fITSResponse=other.fITSResponse;
 150     fTPCResponse=other.fTPCResponse;
 151     fTRDResponse=other.fTRDResponse;
 152     fTOFResponse=other.fTOFResponse;
 153     fEMCALResponse=other.fEMCALResponse;
 154     fRange=other.fRange;
 155     fITSPIDmethod=other.fITSPIDmethod;
 156     fOADBPath=other.fOADBPath;
 157     fCustomTPCpidResponse=other.fCustomTPCpidResponse;
 158     fIsMC=other.fIsMC;
 159     fBeamType="PP";
 160     fLHCperiod="";
 161     fMCperiodTPC="";
 162     fMCperiodUser=other.fMCperiodUser;
 163     fCurrentFile="";
 164     fRecoPass=0;
 165     fRecoPassUser=other.fRecoPassUser;
 166     fRun=0;
 167     fOldRun=0;
 168     fArrPidResponseMaster=0x0;
 169     fResolutionCorrection=0x0;
 170     fTRDPIDResponseObject=0x0;
 171     fEMCALPIDParams=0x0;
 172     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 173     fTOFtail=1.1;
 174     fTOFPIDParams=0x0;
 175     fCurrentEvent=other.fCurrentEvent;
 176   }
 177   return *this;
 178 }
 179
 180 //______________________________________________________________________________
 181 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
 182 {
 183   //
 184   // NumberOfSigmas for 'detCode'
 185   //
 186
 187   switch (detCode){
 188     case kDetITS: return NumberOfSigmasITS(track, type); break;
 189     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
 190     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
 191     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
 192     default: return -999.;
 193   }
 194
 195 }
 196
 197 //______________________________________________________________________________
 198 Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
 199 {
 200   //
 201   // NumberOfSigmas for 'detCode'
 202   //
 203   return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
 204 }
 205
 206 //______________________________________________________________________________
 207 Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
 208 {
 209   //
 210   // Calculate the number of sigmas in the ITS
 211   //
 212
 213   AliVTrack *track=(AliVTrack*)vtrack;
 214
 215   // look for cached value first
 216   // only the non SA tracks are cached
 217   if ( track->GetDetectorPID() ){
 218     return track->GetDetectorPID()->GetNumberOfSigmas(kITS, type);
 219   }
 220
 221   Float_t dEdx=track->GetITSsignal();
 222   if (dEdx<=0) return -999.;
 223
 224   UChar_t clumap=track->GetITSClusterMap();
 225   Int_t nPointsForPid=0;
 226   for(Int_t i=2; i<6; i++){
 227     if(clumap&(1<<i)) ++nPointsForPid;
 228   }
 229   Float_t mom=track->P();
 230
 231   //check for ITS standalone tracks
 232   Bool_t isSA=kTRUE;
 233   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 234
 235   //TODO: in case of the electron, use the SA parametrisation,
 236   //      this needs to be changed if ITS provides a parametrisation
 237   //      for electrons also for ITS+TPC tracks
 238   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
 239 }
 240
 241 //______________________________________________________________________________
 242 Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
 243 {
 244   //
 245   // Calculate the number of sigmas in the TPC
 246   //
 247
 248   AliVTrack *track=(AliVTrack*)vtrack;
 249
 250   // look for cached value first
 251   if (track->GetDetectorPID()){
 252     return track->GetDetectorPID()->GetNumberOfSigmas(kTPC, type);
 253   }
 254
 255   Double_t mom  = track->GetTPCmomentum();
 256   Double_t sig  = track->GetTPCsignal();
 257   UInt_t   sigN = track->GetTPCsignalN();
 258
 259   Double_t nSigma = -999.;
 260   if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
 261
 262   return nSigma;
 263 }
 264
 265 //______________________________________________________________________________
 266 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
 267 {
 268   //
 269   // Calculate the number of sigmas in the EMCAL
 270   //
 271
 272   AliVTrack *track=(AliVTrack*)vtrack;
 273
 274   // look for cached value first
 275   if (track->GetDetectorPID()){
 276     return track->GetDetectorPID()->GetNumberOfSigmas(kEMCAL, type);
 277   }
 278
 279   AliVCluster *matchedClus = NULL;
 280
 281   Double_t mom     = -1.;
 282   Double_t pt      = -1.;
 283   Double_t EovP    = -1.;
 284   Double_t fClsE   = -1.;
 285
 286   Int_t nMatchClus = -1;
 287   Int_t charge     = 0;
 288
 289   // Track matching
 290   nMatchClus = track->GetEMCALcluster();
 291   if(nMatchClus > -1){
 292
 293     mom    = track->P();
 294     pt     = track->Pt();
 295     charge = track->Charge();
 296
 297     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 298
 299     if(matchedClus){
 300
 301       // matched cluster is EMCAL
 302       if(matchedClus->IsEMCAL()){
 303
 304         fClsE       = matchedClus->E();
 305         EovP        = fClsE/mom;
 306
 307
 308         // NSigma value really meaningful only for electrons!
 309         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
 310       }
 311     }
 312   }
 313
 314   return -999;
 315
 316 }
 317
 318 //______________________________________________________________________________
 319 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
 320
 321   AliVTrack *track=(AliVTrack*)vtrack;
 322
 323   AliVCluster *matchedClus = NULL;
 324
 325   Double_t mom     = -1.;
 326   Double_t pt      = -1.;
 327   Double_t EovP    = -1.;
 328   Double_t fClsE   = -1.;
 329
 330   // initialize eop and shower shape parameters
 331   eop = -1.;
 332   for(Int_t i = 0; i < 4; i++){
 333     showershape[i] = -1.;
 334   }
 335
 336   Int_t nMatchClus = -1;
 337   Int_t charge     = 0;
 338
 339   // Track matching
 340   nMatchClus = track->GetEMCALcluster();
 341   if(nMatchClus > -1){
 342
 343     mom    = track->P();
 344     pt     = track->Pt();
 345     charge = track->Charge();
 346
 347     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 348
 349     if(matchedClus){
 350
 351       // matched cluster is EMCAL
 352       if(matchedClus->IsEMCAL()){
 353
 354         fClsE       = matchedClus->E();
 355         EovP        = fClsE/mom;
 356
 357         // fill used EMCAL variables here
 358         eop            = EovP; // E/p
 359         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
 360         showershape[1] = matchedClus->GetM02(); // long axis
 361         showershape[2] = matchedClus->GetM20(); // short axis
 362         showershape[3] = matchedClus->GetDispersion(); // dispersion
 363
 364         // NSigma value really meaningful only for electrons!
 365         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
 366       }
 367     }
 368   }
 369   return -999;
 370 }
 371
 372
 373 //______________________________________________________________________________
 374 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 375 {
 376   //
 377   // Compute PID response of 'detCode'
 378   //
 379
 380   switch (detCode){
 381     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
 382     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
 383     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
 384     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
 385     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
 386     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
 387     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
 388     default: return kDetNoSignal;
 389   }
 390 }
 391
 392 //______________________________________________________________________________
 393 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 394 {
 395   //
 396   // Compute PID response of 'detCode'
 397   //
 398
 399   return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
 400 }
 401
 402 //______________________________________________________________________________
 403 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 404 {
 405   //
 406   // Compute PID response for the ITS
 407   //
 408
 409   // look for cached value first
 410   // only the non SA tracks are cached
 411   if (track->GetDetectorPID()){
 412     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
 413   }
 414
 415   // set flat distribution (no decision)
 416   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 417
 418   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
 419     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
 420
 421   //check for ITS standalone tracks
 422   Bool_t isSA=kTRUE;
 423   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 424
 425   Double_t mom=track->P();
 426   Double_t dedx=track->GetITSsignal();
 427   Double_t momITS=mom;
 428   UChar_t clumap=track->GetITSClusterMap();
 429   Int_t nPointsForPid=0;
 430   for(Int_t i=2; i<6; i++){
 431     if(clumap&(1<<i)) ++nPointsForPid;
 432   }
 433
 434   if(nPointsForPid<3) { // track not to be used for combined PID purposes
 435     //       track->ResetStatus(AliVTrack::kITSpid);
 436     return kDetNoSignal;
 437   }
 438
 439   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
 440   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 441     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
 442     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
 443     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
 444     //TODO: in case of the electron, use the SA parametrisation,
 445     //      this needs to be changed if ITS provides a parametrisation
 446     //      for electrons also for ITS+TPC tracks
 447     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
 448     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 449       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 450     } else {
 451       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 452       mismatch=kFALSE;
 453     }
 454
 455     // Check for particles heavier than (AliPID::kSPECIES - 1)
 456     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
 457
 458   }
 459
 460   if (mismatch){
 461     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 462     return kDetNoSignal;
 463   }
 464
 465
 466   return kDetPidOk;
 467 }
 468 //______________________________________________________________________________
 469 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 470 {
 471   //
 472   // Compute PID response for the TPC
 473   //
 474
 475   // look for cached value first
 476   if (track->GetDetectorPID()){
 477     return track->GetDetectorPID()->GetRawProbability(kTPC, p, nSpecies);
 478   }
 479
 480   // set flat distribution (no decision)
 481   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 482
 483   // check quality of the track
 484   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
 485
 486   Double_t mom = track->GetTPCmomentum();
 487
 488   Double_t dedx=track->GetTPCsignal();
 489   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
 490
 491   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
 492
 493   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
 494     AliPID::EParticleType type=AliPID::EParticleType(j);
 495     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 496     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
 497     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 498       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 499     } else {
 500       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 501       mismatch=kFALSE;
 502     }
 503   }
 504
 505   if (mismatch){
 506     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 507     return kDetNoSignal;
 508   }
 509
 510   return kDetPidOk;
 511 }
 512 //______________________________________________________________________________
 513 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 514 {
 515   //
 516   // Compute PID response for the
 517   //
 518
 519   // look for cached value first
 520   if (track->GetDetectorPID()){
 521     return track->GetDetectorPID()->GetRawProbability(kTOF, p, nSpecies);
 522   }
 523
 524   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
 525
 526   // set flat distribution (no decision)
 527   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 528
 529   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
 530   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
 531
 532   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
 533   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
 534     AliPID::EParticleType type=AliPID::EParticleType(j);
 535     Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
 536
 537     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
 538     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
 539     if (TMath::Abs(nsigmas) > (fRange+2)) {
 540       if(nsigmas < fTOFtail)
 541         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
 542       else
 543         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
 544     } else{
 545       if(nsigmas < fTOFtail)
 546         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
 547       else
 548         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
 549     }
 550
 551     if (TMath::Abs(nsigmas)<5.){
 552       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
 553       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
 554     }
 555   }
 556
 557   if (mismatch){
 558     return kDetMismatch;
 559   }
 560
 561     // TODO: Light nuclei
 562
 563   return kDetPidOk;
 564 }
 565 //______________________________________________________________________________
 566 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 567 {
 568   //
 569   // Compute PID response for the
 570   //
 571
 572   // look for cached value first
 573   if (track->GetDetectorPID()){
 574     return track->GetDetectorPID()->GetRawProbability(kTRD, p, nSpecies);
 575   }
 576
 577   // set flat distribution (no decision)
 578   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 579   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
 580
 581   Float_t mom[6]={0.};
 582   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
 583   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
 584   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
 585   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
 586     mom[ilayer] = track->GetTRDmomentum(ilayer);
 587     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
 588       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
 589     }
 590   }
 591   fTRDResponse.GetResponse(nslices, dedx, mom, p);
 592   return kDetPidOk;
 593 }
 594 //______________________________________________________________________________
 595 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 596 {
 597   //
 598   // Compute PID response for the EMCAL
 599   //
 600
 601   // look for cached value first
 602   if (track->GetDetectorPID()){
 603     return track->GetDetectorPID()->GetRawProbability(kEMCAL, p, nSpecies);
 604   }
 605
 606   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 607
 608   AliVCluster *matchedClus = NULL;
 609
 610   Double_t mom     = -1.;
 611   Double_t pt      = -1.;
 612   Double_t EovP    = -1.;
 613   Double_t fClsE   = -1.;
 614
 615   Int_t nMatchClus = -1;
 616   Int_t charge     = 0;
 617
 618   // Track matching
 619   nMatchClus = track->GetEMCALcluster();
 620
 621   if(nMatchClus > -1){
 622
 623     mom    = track->P();
 624     pt     = track->Pt();
 625     charge = track->Charge();
 626
 627     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 628
 629     if(matchedClus){
 630
 631       // matched cluster is EMCAL
 632       if(matchedClus->IsEMCAL()){
 633
 634       fClsE       = matchedClus->E();
 635       EovP        = fClsE/mom;
 636
 637
 638       // compute the probabilities
 639         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
 640
 641           // in case everything is OK
 642               return kDetPidOk;
 643         }
 644       }
 645     }
 646   }
 647
 648   // in all other cases set flat distribution (no decision)
 649   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
 650   return kDetNoSignal;
 651
 652 }
 653 //______________________________________________________________________________
 654 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
 655 {
 656   //
 657   // Compute PID response for the PHOS
 658   //
 659
 660   // look for cached value first
 661 //   if (track->GetDetectorPID()){
 662 //     return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
 663 //   }
 664
 665   // set flat distribution (no decision)
 666   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 667   return kDetNoSignal;
 668 }
 669 //______________________________________________________________________________
 670 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 671 {
 672   //
 673   // Compute PID response for the HMPID
 674   //
 675
 676
 677   // look for cached value first
 678   if (track->GetDetectorPID()){
 679     return track->GetDetectorPID()->GetRawProbability(kHMPID, p, nSpecies);
 680   }
 681
 682   // set flat distribution (no decision)
 683   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 684   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
 685
 686   track->GetHMPIDpid(p);
 687
 688   return kDetPidOk;
 689 }
 690
 691 //______________________________________________________________________________
 692 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
 693 {
 694   //
 695   // Apply settings for the current event
 696   //
 697   fRecoPass=pass;
 698
 699   fCurrentEvent=0x0;
 700   if (!event) return;
 701   fCurrentEvent=event;
 702   if (run>0) fRun=run;
 703   else fRun=event->GetRunNumber();
 704
 705   if (fRun!=fOldRun){
 706     ExecNewRun();
 707     fOldRun=fRun;
 708   }
 709
 710   //TPC resolution parametrisation PbPb
 711   if ( fResolutionCorrection ){
 712     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
 713     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
 714   }
 715
 716   //TOF resolution
 717   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 718
 719
 720   // Get and set centrality
 721   AliCentrality *centrality = event->GetCentrality();
 722   if(centrality){
 723     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
 724   }
 725   else{
 726     fCurrCentrality = -1;
 727   }
 728 }
 729
 730 //______________________________________________________________________________
 731 void AliPIDResponse::ExecNewRun()
 732 {
 733   //
 734   // Things to Execute upon a new run
 735   //
 736   SetRecoInfo();
 737
 738   SetITSParametrisation();
 739
 740   SetTPCPidResponseMaster();
 741   SetTPCParametrisation();
 742
 743   SetTRDPidResponseMaster();
 744   InitializeTRDResponse();
 745
 746   SetEMCALPidResponseMaster();
 747   InitializeEMCALResponse();
 748
 749   SetTOFPidResponseMaster();
 750   InitializeTOFResponse();
 751 }
 752
 753 //_____________________________________________________
 754 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
 755 {
 756   //
 757   // Get TPC multiplicity in bins of 150
 758   //
 759
 760   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
 761   Double_t tpcMulti=0.;
 762   if(vertexTPC){
 763     Double_t vertexContribTPC=vertexTPC->GetNContributors();
 764     tpcMulti=vertexContribTPC/150.;
 765     if (tpcMulti>20.) tpcMulti=20.;
 766   }
 767
 768   return tpcMulti;
 769 }
 770
 771 //______________________________________________________________________________
 772 void AliPIDResponse::SetRecoInfo()
 773 {
 774   //
 775   // Set reconstruction information
 776   //
 777
 778   //reset information
 779   fLHCperiod="";
 780   fMCperiodTPC="";
 781
 782   fBeamType="";
 783
 784   fBeamType="PP";
 785
 786
 787   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
 788   TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
 789
 790   //find the period by run number (UGLY, but not stored in ESD and AOD... )
 791   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
 792   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
 793   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
 794   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 795   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 796   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 797   else if (fRun>=136851&&fRun<=139517) {
 798     fLHCperiod="LHC10H";
 799     fMCperiodTPC="LHC10H8";
 800     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 801     fBeamType="PBPB";
 802   }
 803   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 804   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
 805   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 806   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 807   // also for 11e,f use 11d
 808   else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 809   else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 810
 811   else if (fRun>=166529) {
 812     fLHCperiod="LHC11H";
 813     fMCperiodTPC="LHC11A10";
 814     fBeamType="PBPB";
 815     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
 816   }
 817
 818
 819   //exception new pp MC productions from 2011
 820   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 821   // exception for 11f1
 822   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 823 }
 824
 825 //______________________________________________________________________________
 826 void AliPIDResponse::SetITSParametrisation()
 827 {
 828   //
 829   // Set the ITS parametrisation
 830   //
 831 }
 832
 833 //______________________________________________________________________________
 834 void AliPIDResponse::SetTPCPidResponseMaster()
 835 {
 836   //
 837   // Load the TPC pid response functions from the OADB
 838   //
 839   //don't load twice for the moment
 840    if (fArrPidResponseMaster) return;
 841
 842
 843   //reset the PID response functions
 844   delete fArrPidResponseMaster;
 845   fArrPidResponseMaster=0x0;
 846
 847   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
 848   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
 849
 850   TFile *f=TFile::Open(fileName.Data());
 851   if (f && f->IsOpen() && !f->IsZombie()){
 852     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
 853   }
 854   delete f;
 855
 856   if (!fArrPidResponseMaster){
 857     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
 858     return;
 859   }
 860   fArrPidResponseMaster->SetOwner();
 861 }
 862
 863 //______________________________________________________________________________
 864 void AliPIDResponse::SetTPCParametrisation()
 865 {
 866   //
 867   // Change BB parametrisation for current run
 868   //
 869
 870   if (fLHCperiod.IsNull()) {
 871     AliFatal("No period set, not changing parametrisation");
 872     return;
 873   }
 874
 875   //
 876   // Set default parametrisations for data and MC
 877   //
 878
 879   //data type
 880   TString datatype="DATA";
 881   //in case of mc fRecoPass is per default 1
 882   if (fIsMC) {
 883       if(!fTuneMConData) datatype="MC";
 884       fRecoPass=1;
 885   }
 886
 887   //
 888   //reset old splines
 889   //
 890   for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
 891     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
 892   }
 893
 894   // period
 895   TString period=fLHCperiod;
 896   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
 897
 898   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 899   Bool_t found=kFALSE;
 900   //
 901   //set the new PID splines
 902   //
 903   if (fArrPidResponseMaster){
 904     Int_t recopass = fRecoPass;
 905     if(fTuneMConData) recopass = fRecoPassUser;
 906     TObject *grAll=0x0;
 907     //for MC don't use period information
 908 //     if (fIsMC) period="[A-Z0-9]*";
 909     //for MC use MC period information
 910 //pattern for the default entry (valid for all particles)
 911     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
 912
 913     //loop over entries and filter them
 914     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
 915       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
 916       if (responseFunction==0x0) continue;
 917       TString responseName=responseFunction->GetName();
 918
 919       if (!reg.MatchB(responseName)) continue;
 920
 921       TObjArray *arr=reg.MatchS(responseName);
 922       TString particleName=arr->At(1)->GetName();
 923       delete arr;
 924       if (particleName.IsNull()) continue;
 925       if (!grAll&&particleName=="ALL") grAll=responseFunction;
 926       else {
 927         //find particle id
 928         for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
 929           TString particle=AliPID::ParticleName(ispec);
 930           particle.ToUpper();
 931           if ( particle == particleName ){
 932             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
 933             fTPCResponse.SetUseDatabase(kTRUE);
 934             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
 935             // overwrite default with proton spline (for light nuclei)
 936             if (ispec==AliPID::kProton) grAll=responseFunction;
 937             found=kTRUE;
 938             break;
 939           }
 940         }
 941       }
 942     }
 943
 944     //set default response function to all particles which don't have a specific one
 945     if (grAll){
 946       for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
 947         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
 948           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
 949           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
 950           found=kTRUE;
 951         }
 952       }
 953     }
 954   }
 955
 956   if (!found){
 957     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 958   }
 959   //
 960   // Setup resolution parametrisation
 961   //
 962
 963   //default
 964   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
 965
 966   if (fRun>=122195){
 967     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
 968   }
 969   if (fArrPidResponseMaster)
 970   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 971
 972   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
 973 }
 974
 975 //______________________________________________________________________________
 976 void AliPIDResponse::SetTRDPidResponseMaster()
 977 {
 978   //
 979   // Load the TRD pid params and references from the OADB
 980   //
 981   if(fTRDPIDResponseObject) return;
 982   AliOADBContainer contParams("contParams");
 983
 984   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
 985   if(statusResponse){
 986     AliError("Failed initializing PID Response Object from OADB");
 987   } else {
 988     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
 989     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
 990     if(!fTRDPIDResponseObject){
 991       AliError(Form("TRD Response not found in run %d", fRun));
 992     }
 993   }
 994   /*
 995   AliOADBContainer contRefs("contRefs");
 996   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
 997   if(statusRefs){
 998     AliInfo("Failed Loading References for TRD");
 999   } else {
1000     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
1001     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
1002     if(!fTRDPIDReference){
1003       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
1004     }
1005     }
1006     */
1007 }
1008
1009 //______________________________________________________________________________
1010 void AliPIDResponse::InitializeTRDResponse(){
1011   //
1012   // Set PID Params and references to the TRD PID response
1013   //
1014     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1015     SetTRDPIDmethod();
1016 }
1017
1018 void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
1019
1020   fTRDResponse.SetPIDmethod(method);
1021   if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1022     // backward compatibility for setting with 8 slices
1023     fTRDslicesForPID[0] = 0;
1024     fTRDslicesForPID[1] = 7;
1025   }
1026   else{
1027     if(method==AliTRDPIDResponse::kLQ1D){
1028       fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1029       fTRDslicesForPID[1] = 0;
1030     }
1031     if(method==AliTRDPIDResponse::kLQ2D){
1032       fTRDslicesForPID[0] = 1;
1033       fTRDslicesForPID[1] = 7;
1034     }
1035   }
1036   AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
1037 }
1038
1039 //______________________________________________________________________________
1040 void AliPIDResponse::SetTOFPidResponseMaster()
1041 {
1042   //
1043   // Load the TOF pid params from the OADB
1044   //
1045
1046   if (fTOFPIDParams) delete fTOFPIDParams;
1047   fTOFPIDParams=0x0;
1048
1049   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1050   if (oadbf && oadbf->IsOpen()) {
1051     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1052     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1053     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1054     oadbf->Close();
1055     delete oadbc;
1056   }
1057   delete oadbf;
1058
1059   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1060 }
1061
1062 //______________________________________________________________________________
1063 void AliPIDResponse::InitializeTOFResponse(){
1064   //
1065   // Set PID Params to the TOF PID response
1066   //
1067
1068   AliInfo("TOF PID Params loaded from OADB");
1069   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1070   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1071   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1072                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1073
1074   for (Int_t i=0;i<4;i++) {
1075     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1076   }
1077   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1078
1079 }
1080
1081
1082 //_________________________________________________________________________
1083 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
1084   //
1085   // Check whether track is identified as electron under a given electron efficiency hypothesis
1086   //
1087   Double_t probs[AliPID::kSPECIES];
1088   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
1089
1090   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1091   // Take mean of the TRD momenta in the given tracklets
1092   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1093   Int_t nmomenta = 0;
1094   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1095     if(vtrack->GetTRDmomentum(iPl) > 0.){
1096       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1097     }
1098   }
1099   p = TMath::Mean(nmomenta, trdmomenta);
1100
1101   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
1102 }
1103
1104 //______________________________________________________________________________
1105 void AliPIDResponse::SetEMCALPidResponseMaster()
1106 {
1107   //
1108   // Load the EMCAL pid response functions from the OADB
1109   //
1110   TObjArray* fEMCALPIDParamsRun      = NULL;
1111   TObjArray* fEMCALPIDParamsPass     = NULL;
1112
1113   if(fEMCALPIDParams) return;
1114   AliOADBContainer contParams("contParams");
1115
1116   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1117   if(statusPars){
1118     AliError("Failed initializing PID Params from OADB");
1119   }
1120   else {
1121     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1122
1123     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1124     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1125     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1126
1127     if(!fEMCALPIDParams){
1128       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1129       AliInfo("Will take the standard LHC11d instead ...");
1130
1131       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1132       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1133       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1134
1135       if(!fEMCALPIDParams){
1136         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1137       }
1138     }
1139   }
1140 }
1141
1142 //______________________________________________________________________________
1143 void AliPIDResponse::InitializeEMCALResponse(){
1144   //
1145   // Set PID Params to the EMCAL PID response
1146   //
1147   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1148
1149 }
1150
1151 //_____________________________________________________
1152 void AliPIDResponse::FillTrackDetectorPID()
1153 {
1154   //
1155   // create detector PID information and setup the transient pointer in the track
1156   //
1157
1158   if (!fCurrentEvent) return;
1159
1160   //TODO: which particles to include? See also the loops below...
1161   Double_t values[AliPID::kSPECIESC]={0};
1162
1163   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1164     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1165     if (!track) continue;
1166
1167     AliDetectorPID *detPID=new AliDetectorPID;
1168     for (Int_t idet=0; idet<kNdetectors; ++idet){
1169
1170       //nsigmas
1171       for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1172         values[ipart]=NumberOfSigmas((EDetector)idet,track,(AliPID::EParticleType)ipart);
1173       detPID->SetNumberOfSigmas((EDetector)idet, values, (Int_t)AliPID::kSPECIESC);
1174
1175       //probabilities
1176       EDetPidStatus status=ComputePIDProbability((EDetector)idet,track,AliPID::kSPECIESC,values);
1177       detPID->SetRawProbability((EDetector)idet, values, (Int_t)AliPID::kSPECIESC, status);
1178     }
1179
1180     track->SetDetectorPID(detPID);
1181   }
1182 }
1183
1184 //_________________________________________________________________________
1185 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1186   //
1187   // Set TOF response function
1188   // Input option for event_time used
1189   //
1190
1191     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1192     if(t0spread < 10) t0spread = 80;
1193
1194     // T0 from TOF algorithm
1195
1196     Bool_t flagT0TOF=kFALSE;
1197     Bool_t flagT0T0=kFALSE;
1198     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1199     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1200     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1201
1202     // T0-TOF arrays
1203     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1204     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1205     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1206       estimatedT0event[i]=0.0;
1207       estimatedT0resolution[i]=0.0;
1208       startTimeMask[i] = 0;
1209     }
1210
1211     Float_t resT0A=75,resT0C=65,resT0AC=55;
1212     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1213         flagT0T0=kTRUE;
1214     }
1215
1216
1217     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1218
1219     if (tofHeader) { // read global info and T0-TOF
1220       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1221       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1222       if(t0spread < 10) t0spread = 80;
1223
1224       flagT0TOF=kTRUE;
1225       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1226         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1227         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1228         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1229       }
1230
1231       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1232       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1233       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1234       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1235         Int_t icurrent = (Int_t)ibin->GetAt(j);
1236         startTime[icurrent]=t0Bin->GetAt(j);
1237         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1238         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1239       }
1240     }
1241
1242     // for cut of 3 sigma on t0 spread
1243     Float_t t0cut = 3 * t0spread;
1244     if(t0cut < 500) t0cut = 500;
1245
1246     if(option == kFILL_T0){ // T0-FILL is used
1247         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1248           estimatedT0event[i]=0.0;
1249           estimatedT0resolution[i]=t0spread;
1250         }
1251         fTOFResponse.SetT0event(estimatedT0event);
1252         fTOFResponse.SetT0resolution(estimatedT0resolution);
1253     }
1254
1255     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1256         if(flagT0TOF){
1257             fTOFResponse.SetT0event(startTime);
1258             fTOFResponse.SetT0resolution(startTimeRes);
1259             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1260               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1261               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1262             }
1263         }
1264         else{
1265             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1266               estimatedT0event[i]=0.0;
1267               estimatedT0resolution[i]=t0spread;
1268               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1269             }
1270             fTOFResponse.SetT0event(estimatedT0event);
1271             fTOFResponse.SetT0resolution(estimatedT0resolution);
1272         }
1273     }
1274     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1275         Float_t t0AC=-10000;
1276         Float_t t0A=-10000;
1277         Float_t t0C=-10000;
1278         if(flagT0T0){
1279             t0AC= vevent->GetT0TOF()[0];
1280             t0A= vevent->GetT0TOF()[1];
1281             t0C= vevent->GetT0TOF()[2];
1282         }
1283
1284         Float_t t0t0Best = 0;
1285         Float_t t0t0BestRes = 9999;
1286         Int_t t0used=0;
1287         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1288             t0t0Best = t0AC;
1289             t0t0BestRes = resT0AC;
1290             t0used=6;
1291         }
1292         else if(TMath::Abs(t0C) < t0cut){
1293             t0t0Best = t0C;
1294             t0t0BestRes = resT0C;
1295             t0used=4;
1296         }
1297         else if(TMath::Abs(t0A) < t0cut){
1298             t0t0Best = t0A;
1299             t0t0BestRes = resT0A;
1300             t0used=2;
1301         }
1302
1303         if(flagT0TOF){ // if T0-TOF info is available
1304             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1305                 if(t0t0BestRes < 999){
1306                   if(startTimeRes[i] < t0spread){
1307                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1308                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1309                     estimatedT0event[i]=t0best / wtot;
1310                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1311                     startTimeMask[i] = t0used+1;
1312                   }
1313                   else {
1314                     estimatedT0event[i]=t0t0Best;
1315                     estimatedT0resolution[i]=t0t0BestRes;
1316                     startTimeMask[i] = t0used;
1317                   }
1318                 }
1319                 else{
1320                   estimatedT0event[i]=startTime[i];
1321                   estimatedT0resolution[i]=startTimeRes[i];
1322                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1323                 }
1324                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1325             }
1326             fTOFResponse.SetT0event(estimatedT0event);
1327             fTOFResponse.SetT0resolution(estimatedT0resolution);
1328         }
1329         else{ // if no T0-TOF info is available
1330             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1331               fTOFResponse.SetT0binMask(i,t0used);
1332               if(t0t0BestRes < 999){
1333                 estimatedT0event[i]=t0t0Best;
1334                 estimatedT0resolution[i]=t0t0BestRes;
1335               }
1336               else{
1337                 estimatedT0event[i]=0.0;
1338                 estimatedT0resolution[i]=t0spread;
1339               }
1340             }
1341             fTOFResponse.SetT0event(estimatedT0event);
1342             fTOFResponse.SetT0resolution(estimatedT0resolution);
1343         }
1344     }
1345
1346     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1347         Float_t t0AC=-10000;
1348         Float_t t0A=-10000;
1349         Float_t t0C=-10000;
1350         if(flagT0T0){
1351             t0AC= vevent->GetT0TOF()[0];
1352             t0A= vevent->GetT0TOF()[1];
1353             t0C= vevent->GetT0TOF()[2];
1354         }
1355
1356         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1357             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1358               estimatedT0event[i]=t0AC;
1359               estimatedT0resolution[i]=resT0AC;
1360               fTOFResponse.SetT0binMask(i,6);
1361             }
1362         }
1363         else if(TMath::Abs(t0C) < t0cut){
1364             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1365               estimatedT0event[i]=t0C;
1366               estimatedT0resolution[i]=resT0C;
1367               fTOFResponse.SetT0binMask(i,4);
1368             }
1369         }
1370         else if(TMath::Abs(t0A) < t0cut){
1371             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1372               estimatedT0event[i]=t0A;
1373               estimatedT0resolution[i]=resT0A;
1374               fTOFResponse.SetT0binMask(i,2);
1375             }
1376         }
1377         else{
1378             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1379               estimatedT0event[i]=0.0;
1380               estimatedT0resolution[i]=t0spread;
1381               fTOFResponse.SetT0binMask(i,0);
1382             }
1383         }
1384         fTOFResponse.SetT0event(estimatedT0event);
1385         fTOFResponse.SetT0resolution(estimatedT0resolution);
1386     }
1387     delete [] startTime;
1388     delete [] startTimeRes;
1389     delete [] startTimeMask;
1390     delete [] estimatedT0event;
1391     delete [] estimatedT0resolution;
1392 }