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
 836   //exception new pp MC productions from 2011
 837   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 838   // exception for 11f1
 839   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 840 }
 841
 842 //______________________________________________________________________________
 843 void AliPIDResponse::SetITSParametrisation()
 844 {
 845   //
 846   // Set the ITS parametrisation
 847   //
 848 }
 849
 850 //______________________________________________________________________________
 851 void AliPIDResponse::SetTPCPidResponseMaster()
 852 {
 853   //
 854   // Load the TPC pid response functions from the OADB
 855   // Load the TPC voltage maps from OADB
 856   //
 857   //don't load twice for the moment
 858    if (fArrPidResponseMaster) return;
 859
 860
 861   //reset the PID response functions
 862   delete fArrPidResponseMaster;
 863   fArrPidResponseMaster=NULL;
 864
 865   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
 866   TFile *f=NULL;
 867   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
 868
 869   TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
 870   f=TFile::Open(fileNamePIDresponse.Data());
 871   if (f && f->IsOpen() && !f->IsZombie()){
 872     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
 873   }
 874   delete f;
 875
 876   TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
 877   f=TFile::Open(fileNameVoltageMaps.Data());
 878   if (f && f->IsOpen() && !f->IsZombie()){
 879     fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
 880   }
 881   delete f;
 882
 883   if (!fArrPidResponseMaster){
 884     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
 885     return;
 886   }
 887   fArrPidResponseMaster->SetOwner();
 888
 889   if (!fOADBvoltageMaps)
 890   {
 891     AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
 892   }
 893   fArrPidResponseMaster->SetOwner();
 894 }
 895
 896 //______________________________________________________________________________
 897 void AliPIDResponse::SetTPCParametrisation()
 898 {
 899   //
 900   // Change BB parametrisation for current run
 901   //
 902
 903   if (fLHCperiod.IsNull()) {
 904     AliFatal("No period set, not changing parametrisation");
 905     return;
 906   }
 907
 908   //
 909   // Set default parametrisations for data and MC
 910   //
 911
 912   //data type
 913   TString datatype="DATA";
 914   //in case of mc fRecoPass is per default 1
 915   if (fIsMC) {
 916       if(!fTuneMConData) datatype="MC";
 917       fRecoPass=1;
 918   }
 919
 920   //
 921   //reset old splines
 922   //
 923   fTPCResponse.ResetSplines();
 924
 925   // period
 926   TString period=fLHCperiod;
 927   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
 928
 929   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 930   Bool_t found=kFALSE;
 931   //
 932   //set the new PID splines
 933   //
 934   if (fArrPidResponseMaster){
 935     Int_t recopass = fRecoPass;
 936     if(fTuneMConData) recopass = fRecoPassUser;
 937     //for MC don't use period information
 938     //if (fIsMC) period="[A-Z0-9]*";
 939     //for MC use MC period information
 940     //pattern for the default entry (valid for all particles)
 941     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
 942
 943     //find particle id ang gain scenario
 944     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
 945     {
 946       TObject *grAll=NULL;
 947       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
 948       gainScenario.ToUpper();
 949       //loop over entries and filter them
 950       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
 951       {
 952         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
 953         if (responseFunction==NULL) continue;
 954         TString responseName=responseFunction->GetName();
 955
 956         if (!reg.MatchB(responseName)) continue;
 957
 958         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
 959         TObject* tmp=NULL;
 960         tmp=arr->At(1); if (!tmp) continue;
 961         TString particleName=tmp->GetName();
 962         tmp=arr->At(3); if (!tmp) continue;
 963         TString gainScenarioName=tmp->GetName();
 964         delete arr;
 965         if (particleName.IsNull()) continue;
 966         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
 967         else
 968         {
 969           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
 970           {
 971             TString particle=AliPID::ParticleName(ispec);
 972             particle.ToUpper();
 973             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
 974             if ( particle == particleName && gainScenario == gainScenarioName )
 975             {
 976               fTPCResponse.SetResponseFunction( responseFunction,
 977                                                 (AliPID::EParticleType)ispec,
 978                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
 979               fTPCResponse.SetUseDatabase(kTRUE);
 980               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
 981               found=kTRUE;
 982               // overwrite default with proton spline (for light nuclei)
 983               if (ispec==AliPID::kProton) grAll=responseFunction;
 984               break;
 985             }
 986           }
 987         }
 988       }
 989       if (grAll)
 990       {
 991         for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
 992         {
 993           if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
 994                                                  (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
 995           {
 996               fTPCResponse.SetResponseFunction( grAll,
 997                                                 (AliPID::EParticleType)ispec,
 998                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
 999               fTPCResponse.SetUseDatabase(kTRUE);
1000               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1001               found=kTRUE;
1002           }
1003         }
1004       }
1005     }
1006   }
1007   else AliInfo("no fArrPidResponseMaster");
1008
1009   if (!found){
1010     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1011   }
1012
1013   //
1014   // Setup resolution parametrisation
1015   //
1016
1017   //default
1018   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1019
1020   if (fRun>=122195){
1021     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1022   }
1023   if (fArrPidResponseMaster)
1024   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1025
1026   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1027
1028   //read in the voltage map
1029   TVectorF* gsm = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1030   if (gsm)
1031   {
1032     fTPCResponse.SetVoltageMap(*gsm);
1033     TString vals;
1034     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1035     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1036     AliInfo(vals.Data());
1037     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1038     AliInfo(vals.Data());
1039     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1040     AliInfo(vals.Data());
1041     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1042     AliInfo(vals.Data());
1043   }
1044   else AliInfo("no voltage map, ideal default assumed");
1045 }
1046
1047 //______________________________________________________________________________
1048 void AliPIDResponse::SetTRDPidResponseMaster()
1049 {
1050   //
1051   // Load the TRD pid params and references from the OADB
1052   //
1053   if(fTRDPIDResponseObject) return;
1054   AliOADBContainer contParams("contParams");
1055
1056   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1057   if(statusResponse){
1058     AliError("Failed initializing PID Response Object from OADB");
1059   } else {
1060     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1061     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1062     if(!fTRDPIDResponseObject){
1063       AliError(Form("TRD Response not found in run %d", fRun));
1064     }
1065   }
1066   /*
1067   AliOADBContainer contRefs("contRefs");
1068   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
1069   if(statusRefs){
1070     AliInfo("Failed Loading References for TRD");
1071   } else {
1072     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
1073     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
1074     if(!fTRDPIDReference){
1075       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
1076     }
1077     }
1078     */
1079 }
1080
1081 //______________________________________________________________________________
1082 void AliPIDResponse::InitializeTRDResponse(){
1083   //
1084   // Set PID Params and references to the TRD PID response
1085   //
1086     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1087 }
1088
1089 //______________________________________________________________________________
1090 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1091
1092     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1093         // backward compatibility for setting with 8 slices
1094         TRDslicesForPID[0] = 0;
1095         TRDslicesForPID[1] = 7;
1096     }
1097     else{
1098         if(method==AliTRDPIDResponse::kLQ1D){
1099             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1100             TRDslicesForPID[1] = 0;
1101         }
1102         if(method==AliTRDPIDResponse::kLQ2D){
1103             TRDslicesForPID[0] = 1;
1104             TRDslicesForPID[1] = 7;
1105         }
1106     }
1107     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1108 }
1109
1110 //______________________________________________________________________________
1111 void AliPIDResponse::SetTOFPidResponseMaster()
1112 {
1113   //
1114   // Load the TOF pid params from the OADB
1115   //
1116
1117   if (fTOFPIDParams) delete fTOFPIDParams;
1118   fTOFPIDParams=NULL;
1119
1120   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1121   if (oadbf && oadbf->IsOpen()) {
1122     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1123     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1124     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1125     oadbf->Close();
1126     delete oadbc;
1127   }
1128   delete oadbf;
1129
1130   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1131 }
1132
1133 //______________________________________________________________________________
1134 void AliPIDResponse::InitializeTOFResponse(){
1135   //
1136   // Set PID Params to the TOF PID response
1137   //
1138
1139   AliInfo("TOF PID Params loaded from OADB");
1140   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1141   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1142   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1143                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1144
1145   for (Int_t i=0;i<4;i++) {
1146     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1147   }
1148   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1149
1150 }
1151
1152
1153 //_________________________________________________________________________
1154 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1155   //
1156   // Check whether track is identified as electron under a given electron efficiency hypothesis
1157     //
1158
1159   Double_t probs[AliPID::kSPECIES];
1160   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1161
1162   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1163   // Take mean of the TRD momenta in the given tracklets
1164   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1165   Int_t nmomenta = 0;
1166   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1167     if(vtrack->GetTRDmomentum(iPl) > 0.){
1168       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1169     }
1170   }
1171   p = TMath::Mean(nmomenta, trdmomenta);
1172
1173   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1174 }
1175
1176 //______________________________________________________________________________
1177 void AliPIDResponse::SetEMCALPidResponseMaster()
1178 {
1179   //
1180   // Load the EMCAL pid response functions from the OADB
1181   //
1182   TObjArray* fEMCALPIDParamsRun      = NULL;
1183   TObjArray* fEMCALPIDParamsPass     = NULL;
1184
1185   if(fEMCALPIDParams) return;
1186   AliOADBContainer contParams("contParams");
1187
1188   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1189   if(statusPars){
1190     AliError("Failed initializing PID Params from OADB");
1191   }
1192   else {
1193     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1194
1195     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1196     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1197     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1198
1199     if(!fEMCALPIDParams){
1200       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1201       AliInfo("Will take the standard LHC11d instead ...");
1202
1203       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1204       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1205       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1206
1207       if(!fEMCALPIDParams){
1208         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1209       }
1210     }
1211   }
1212 }
1213
1214 //______________________________________________________________________________
1215 void AliPIDResponse::InitializeEMCALResponse(){
1216   //
1217   // Set PID Params to the EMCAL PID response
1218   //
1219   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1220
1221 }
1222
1223 //_____________________________________________________
1224 void AliPIDResponse::FillTrackDetectorPID()
1225 {
1226   //
1227   // create detector PID information and setup the transient pointer in the track
1228   //
1229
1230   if (!fCurrentEvent) return;
1231
1232   //TODO: which particles to include? See also the loops below...
1233   Double_t values[AliPID::kSPECIESC]={0};
1234
1235   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1236     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1237     if (!track) continue;
1238
1239     AliDetectorPID *detPID=new AliDetectorPID;
1240     for (Int_t idet=0; idet<kNdetectors; ++idet){
1241
1242       //nsigmas
1243       for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1244         values[ipart]=NumberOfSigmas((EDetector)idet,track,(AliPID::EParticleType)ipart);
1245       detPID->SetNumberOfSigmas((EDetector)idet, values, (Int_t)AliPID::kSPECIESC);
1246
1247       //probabilities
1248       EDetPidStatus status=ComputePIDProbability((EDetector)idet,track,AliPID::kSPECIESC,values);
1249       detPID->SetRawProbability((EDetector)idet, values, (Int_t)AliPID::kSPECIESC, status);
1250     }
1251
1252     track->SetDetectorPID(detPID);
1253   }
1254 }
1255
1256 //_________________________________________________________________________
1257 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1258   //
1259   // Set TOF response function
1260   // Input option for event_time used
1261   //
1262
1263     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1264     if(t0spread < 10) t0spread = 80;
1265
1266     // T0 from TOF algorithm
1267
1268     Bool_t flagT0TOF=kFALSE;
1269     Bool_t flagT0T0=kFALSE;
1270     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1271     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1272     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1273
1274     // T0-TOF arrays
1275     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1276     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1277     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1278       estimatedT0event[i]=0.0;
1279       estimatedT0resolution[i]=0.0;
1280       startTimeMask[i] = 0;
1281     }
1282
1283     Float_t resT0A=75,resT0C=65,resT0AC=55;
1284     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1285         flagT0T0=kTRUE;
1286     }
1287
1288
1289     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1290
1291     if (tofHeader) { // read global info and T0-TOF
1292       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1293       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1294       if(t0spread < 10) t0spread = 80;
1295
1296       flagT0TOF=kTRUE;
1297       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1298         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1299         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1300         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1301       }
1302
1303       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1304       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1305       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1306       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1307         Int_t icurrent = (Int_t)ibin->GetAt(j);
1308         startTime[icurrent]=t0Bin->GetAt(j);
1309         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1310         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1311       }
1312     }
1313
1314     // for cut of 3 sigma on t0 spread
1315     Float_t t0cut = 3 * t0spread;
1316     if(t0cut < 500) t0cut = 500;
1317
1318     if(option == kFILL_T0){ // T0-FILL is used
1319         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1320           estimatedT0event[i]=0.0;
1321           estimatedT0resolution[i]=t0spread;
1322         }
1323         fTOFResponse.SetT0event(estimatedT0event);
1324         fTOFResponse.SetT0resolution(estimatedT0resolution);
1325     }
1326
1327     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1328         if(flagT0TOF){
1329             fTOFResponse.SetT0event(startTime);
1330             fTOFResponse.SetT0resolution(startTimeRes);
1331             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1332               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1333               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1334             }
1335         }
1336         else{
1337             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1338               estimatedT0event[i]=0.0;
1339               estimatedT0resolution[i]=t0spread;
1340               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1341             }
1342             fTOFResponse.SetT0event(estimatedT0event);
1343             fTOFResponse.SetT0resolution(estimatedT0resolution);
1344         }
1345     }
1346     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
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         Float_t t0t0Best = 0;
1357         Float_t t0t0BestRes = 9999;
1358         Int_t t0used=0;
1359         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1360             t0t0Best = t0AC;
1361             t0t0BestRes = resT0AC;
1362             t0used=6;
1363         }
1364         else if(TMath::Abs(t0C) < t0cut){
1365             t0t0Best = t0C;
1366             t0t0BestRes = resT0C;
1367             t0used=4;
1368         }
1369         else if(TMath::Abs(t0A) < t0cut){
1370             t0t0Best = t0A;
1371             t0t0BestRes = resT0A;
1372             t0used=2;
1373         }
1374
1375         if(flagT0TOF){ // if T0-TOF info is available
1376             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1377                 if(t0t0BestRes < 999){
1378                   if(startTimeRes[i] < t0spread){
1379                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1380                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1381                     estimatedT0event[i]=t0best / wtot;
1382                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1383                     startTimeMask[i] = t0used+1;
1384                   }
1385                   else {
1386                     estimatedT0event[i]=t0t0Best;
1387                     estimatedT0resolution[i]=t0t0BestRes;
1388                     startTimeMask[i] = t0used;
1389                   }
1390                 }
1391                 else{
1392                   estimatedT0event[i]=startTime[i];
1393                   estimatedT0resolution[i]=startTimeRes[i];
1394                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1395                 }
1396                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1397             }
1398             fTOFResponse.SetT0event(estimatedT0event);
1399             fTOFResponse.SetT0resolution(estimatedT0resolution);
1400         }
1401         else{ // if no T0-TOF info is available
1402             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1403               fTOFResponse.SetT0binMask(i,t0used);
1404               if(t0t0BestRes < 999){
1405                 estimatedT0event[i]=t0t0Best;
1406                 estimatedT0resolution[i]=t0t0BestRes;
1407               }
1408               else{
1409                 estimatedT0event[i]=0.0;
1410                 estimatedT0resolution[i]=t0spread;
1411               }
1412             }
1413             fTOFResponse.SetT0event(estimatedT0event);
1414             fTOFResponse.SetT0resolution(estimatedT0resolution);
1415         }
1416     }
1417
1418     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1419         Float_t t0AC=-10000;
1420         Float_t t0A=-10000;
1421         Float_t t0C=-10000;
1422         if(flagT0T0){
1423             t0AC= vevent->GetT0TOF()[0];
1424             t0A= vevent->GetT0TOF()[1];
1425             t0C= vevent->GetT0TOF()[2];
1426         }
1427
1428         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1429             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1430               estimatedT0event[i]=t0AC;
1431               estimatedT0resolution[i]=resT0AC;
1432               fTOFResponse.SetT0binMask(i,6);
1433             }
1434         }
1435         else if(TMath::Abs(t0C) < t0cut){
1436             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1437               estimatedT0event[i]=t0C;
1438               estimatedT0resolution[i]=resT0C;
1439               fTOFResponse.SetT0binMask(i,4);
1440             }
1441         }
1442         else if(TMath::Abs(t0A) < t0cut){
1443             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1444               estimatedT0event[i]=t0A;
1445               estimatedT0resolution[i]=resT0A;
1446               fTOFResponse.SetT0binMask(i,2);
1447             }
1448         }
1449         else{
1450             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1451               estimatedT0event[i]=0.0;
1452               estimatedT0resolution[i]=t0spread;
1453               fTOFResponse.SetT0binMask(i,0);
1454             }
1455         }
1456         fTOFResponse.SetT0event(estimatedT0event);
1457         fTOFResponse.SetT0resolution(estimatedT0resolution);
1458     }
1459     delete [] startTime;
1460     delete [] startTimeRes;
1461     delete [] startTimeMask;
1462     delete [] estimatedT0event;
1463     delete [] estimatedT0resolution;
1464 }