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