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