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