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<=139517) {
 682     fLHCperiod="LHC10H";
 683     fMCperiodTPC="LHC10H8";
 684     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 685     fBeamType="PBPB";
 686   }
 687   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 688   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
 689   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 690   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 691   // also for 11e,f use 11d
 692   else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 693   else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 694
 695   else if (fRun>=166529 && 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   if (fLHCperiod.IsNull()) {
 782     AliFatal("No period set, not changing parametrisation");
 783     return;
 784   }
 785
 786   //
 787   // Set default parametrisations for data and MC
 788   //
 789
 790   //data type
 791   TString datatype="DATA";
 792   //in case of mc fRecoPass is per default 1
 793   if (fIsMC) {
 794       if(!fTuneMConData) datatype="MC";
 795       fRecoPass=1;
 796   }
 797
 798   //
 799   //reset old splines
 800   //
 801   fTPCResponse.ResetSplines();
 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 = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
 914   if (gsm)
 915   {
 916     fTPCResponse.SetVoltageMap(*gsm);
 917     TString vals;
 918     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
 919     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
 920     AliInfo(vals.Data());
 921     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
 922     AliInfo(vals.Data());
 923     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
 924     AliInfo(vals.Data());
 925     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
 926     AliInfo(vals.Data());
 927   }
 928   else AliInfo("no voltage map, ideal default assumed");
 929 }
 930
 931 //______________________________________________________________________________
 932 void AliPIDResponse::SetTRDPidResponseMaster()
 933 {
 934   //
 935   // Load the TRD pid params and references from the OADB
 936   //
 937   if(fTRDPIDResponseObject) return;
 938   AliOADBContainer contParams("contParams");
 939
 940   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
 941   if(statusResponse){
 942     AliError("Failed initializing PID Response Object from OADB");
 943   } else {
 944     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
 945     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
 946     if(!fTRDPIDResponseObject){
 947       AliError(Form("TRD Response not found in run %d", fRun));
 948     }
 949   }
 950 }
 951
 952 //______________________________________________________________________________
 953 void AliPIDResponse::InitializeTRDResponse(){
 954   //
 955   // Set PID Params and references to the TRD PID response
 956   //
 957     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
 958 }
 959
 960 //______________________________________________________________________________
 961 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
 962
 963     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
 964         // backward compatibility for setting with 8 slices
 965         TRDslicesForPID[0] = 0;
 966         TRDslicesForPID[1] = 7;
 967     }
 968     else{
 969         if(method==AliTRDPIDResponse::kLQ1D){
 970             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
 971             TRDslicesForPID[1] = 0;
 972         }
 973         if(method==AliTRDPIDResponse::kLQ2D){
 974             TRDslicesForPID[0] = 1;
 975             TRDslicesForPID[1] = 7;
 976         }
 977     }
 978     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
 979 }
 980
 981 //______________________________________________________________________________
 982 void AliPIDResponse::SetTOFPidResponseMaster()
 983 {
 984   //
 985   // Load the TOF pid params from the OADB
 986   //
 987
 988   if (fTOFPIDParams) delete fTOFPIDParams;
 989   fTOFPIDParams=NULL;
 990
 991   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
 992   if (oadbf && oadbf->IsOpen()) {
 993     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
 994     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
 995     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
 996     oadbf->Close();
 997     delete oadbc;
 998   }
 999   delete oadbf;
1000
1001   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1002 }
1003
1004 //______________________________________________________________________________
1005 void AliPIDResponse::InitializeTOFResponse(){
1006   //
1007   // Set PID Params to the TOF PID response
1008   //
1009
1010   AliInfo("TOF PID Params loaded from OADB");
1011   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1012   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1013   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1014                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1015
1016   for (Int_t i=0;i<4;i++) {
1017     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1018   }
1019   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1020
1021   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1022   Float_t t0Spread[4];
1023   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1024   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]));
1025   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1026   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1027   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1028     fResT0AC=t0Spread[3];
1029     fResT0A=TMath::Sqrt(a);
1030     fResT0C=TMath::Sqrt(c);
1031   } else {
1032     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1033     fResT0A=75.;
1034     fResT0C=65.;
1035     fResT0AC=55.;
1036   }
1037   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1038
1039 }
1040
1041
1042 //______________________________________________________________________________
1043 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1044   //
1045   // Check whether track is identified as electron under a given electron efficiency hypothesis
1046     //
1047
1048   Double_t probs[AliPID::kSPECIES];
1049   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1050
1051   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1052   // Take mean of the TRD momenta in the given tracklets
1053   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1054   Int_t nmomenta = 0;
1055   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1056     if(vtrack->GetTRDmomentum(iPl) > 0.){
1057       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1058     }
1059   }
1060   p = TMath::Mean(nmomenta, trdmomenta);
1061
1062   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1063 }
1064
1065 //______________________________________________________________________________
1066 void AliPIDResponse::SetEMCALPidResponseMaster()
1067 {
1068   //
1069   // Load the EMCAL pid response functions from the OADB
1070   //
1071   TObjArray* fEMCALPIDParamsRun      = NULL;
1072   TObjArray* fEMCALPIDParamsPass     = NULL;
1073
1074   if(fEMCALPIDParams) return;
1075   AliOADBContainer contParams("contParams");
1076
1077   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1078   if(statusPars){
1079     AliError("Failed initializing PID Params from OADB");
1080   }
1081   else {
1082     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1083
1084     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1085     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1086     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1087
1088     if(!fEMCALPIDParams){
1089       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1090       AliInfo("Will take the standard LHC11d instead ...");
1091
1092       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1093       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1094       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1095
1096       if(!fEMCALPIDParams){
1097         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1098       }
1099     }
1100   }
1101 }
1102
1103 //______________________________________________________________________________
1104 void AliPIDResponse::InitializeEMCALResponse(){
1105   //
1106   // Set PID Params to the EMCAL PID response
1107   //
1108   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1109
1110 }
1111
1112 //______________________________________________________________________________
1113 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1114 {
1115   //
1116   // create detector PID information and setup the transient pointer in the track
1117   //
1118
1119   // check if detector number is inside accepted range
1120   if (detector == kNdetectors) return;
1121
1122   // get detector pid
1123   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1124   if (!detPID) {
1125     detPID=new AliDetectorPID;
1126     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1127   }
1128
1129   //check if values exist
1130   if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return;
1131
1132   //TODO: which particles to include? See also the loops below...
1133   Double_t values[AliPID::kSPECIESC]={0};
1134
1135   //nsigmas
1136   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1137     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1138   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC);
1139
1140   //probabilities
1141   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1142   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1143 }
1144
1145 //______________________________________________________________________________
1146 void AliPIDResponse::FillTrackDetectorPID()
1147 {
1148   //
1149   // create detector PID information and setup the transient pointer in the track
1150   //
1151
1152   if (!fCurrentEvent) return;
1153
1154   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1155     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1156     if (!track) continue;
1157
1158     for (Int_t idet=0; idet<kNdetectors; ++idet){
1159       FillTrackDetectorPID(track, (EDetector)idet);
1160     }
1161   }
1162 }
1163
1164 //______________________________________________________________________________
1165 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1166   //
1167   // Set TOF response function
1168   // Input option for event_time used
1169   //
1170
1171     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1172     if(t0spread < 10) t0spread = 80;
1173
1174     // T0 from TOF algorithm
1175
1176     Bool_t flagT0TOF=kFALSE;
1177     Bool_t flagT0T0=kFALSE;
1178     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1179     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1180     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1181
1182     // T0-TOF arrays
1183     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1184     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1185     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1186       estimatedT0event[i]=0.0;
1187       estimatedT0resolution[i]=0.0;
1188       startTimeMask[i] = 0;
1189     }
1190
1191     Float_t resT0A=fResT0A;
1192     Float_t resT0C=fResT0C;
1193     Float_t resT0AC=fResT0AC;
1194     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1195         flagT0T0=kTRUE;
1196     }
1197
1198
1199     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1200
1201     if (tofHeader) { // read global info and T0-TOF
1202       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1203       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1204       if(t0spread < 10) t0spread = 80;
1205
1206       flagT0TOF=kTRUE;
1207       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1208         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1209         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1210         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1211       }
1212
1213       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1214       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1215       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1216       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1217         Int_t icurrent = (Int_t)ibin->GetAt(j);
1218         startTime[icurrent]=t0Bin->GetAt(j);
1219         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1220         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1221       }
1222     }
1223
1224     // for cut of 3 sigma on t0 spread
1225     Float_t t0cut = 3 * t0spread;
1226     if(t0cut < 500) t0cut = 500;
1227
1228     if(option == kFILL_T0){ // T0-FILL is used
1229         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1230           estimatedT0event[i]=0.0;
1231           estimatedT0resolution[i]=t0spread;
1232         }
1233         fTOFResponse.SetT0event(estimatedT0event);
1234         fTOFResponse.SetT0resolution(estimatedT0resolution);
1235     }
1236
1237     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1238         if(flagT0TOF){
1239             fTOFResponse.SetT0event(startTime);
1240             fTOFResponse.SetT0resolution(startTimeRes);
1241             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1242               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1243               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1244             }
1245         }
1246         else{
1247             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1248               estimatedT0event[i]=0.0;
1249               estimatedT0resolution[i]=t0spread;
1250               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1251             }
1252             fTOFResponse.SetT0event(estimatedT0event);
1253             fTOFResponse.SetT0resolution(estimatedT0resolution);
1254         }
1255     }
1256     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1257         Float_t t0AC=-10000;
1258         Float_t t0A=-10000;
1259         Float_t t0C=-10000;
1260         if(flagT0T0){
1261             t0A= vevent->GetT0TOF()[1];
1262             t0C= vevent->GetT0TOF()[2];
1263             //      t0AC= vevent->GetT0TOF()[0];
1264             t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1265             resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1266             t0AC /= resT0AC*resT0AC;
1267         }
1268
1269         Float_t t0t0Best = 0;
1270         Float_t t0t0BestRes = 9999;
1271         Int_t t0used=0;
1272         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1273             t0t0Best = t0AC;
1274             t0t0BestRes = resT0AC;
1275             t0used=6;
1276         }
1277         else if(TMath::Abs(t0C) < t0cut){
1278             t0t0Best = t0C;
1279             t0t0BestRes = resT0C;
1280             t0used=4;
1281         }
1282         else if(TMath::Abs(t0A) < t0cut){
1283             t0t0Best = t0A;
1284             t0t0BestRes = resT0A;
1285             t0used=2;
1286         }
1287
1288         if(flagT0TOF){ // if T0-TOF info is available
1289             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1290                 if(t0t0BestRes < 999){
1291                   if(startTimeRes[i] < t0spread){
1292                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1293                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1294                     estimatedT0event[i]=t0best / wtot;
1295                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1296                     startTimeMask[i] = t0used+1;
1297                   }
1298                   else {
1299                     estimatedT0event[i]=t0t0Best;
1300                     estimatedT0resolution[i]=t0t0BestRes;
1301                     startTimeMask[i] = t0used;
1302                   }
1303                 }
1304                 else{
1305                   estimatedT0event[i]=startTime[i];
1306                   estimatedT0resolution[i]=startTimeRes[i];
1307                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1308                 }
1309                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1310             }
1311             fTOFResponse.SetT0event(estimatedT0event);
1312             fTOFResponse.SetT0resolution(estimatedT0resolution);
1313         }
1314         else{ // if no T0-TOF info is available
1315             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1316               fTOFResponse.SetT0binMask(i,t0used);
1317               if(t0t0BestRes < 999){
1318                 estimatedT0event[i]=t0t0Best;
1319                 estimatedT0resolution[i]=t0t0BestRes;
1320               }
1321               else{
1322                 estimatedT0event[i]=0.0;
1323                 estimatedT0resolution[i]=t0spread;
1324               }
1325             }
1326             fTOFResponse.SetT0event(estimatedT0event);
1327             fTOFResponse.SetT0resolution(estimatedT0resolution);
1328         }
1329     }
1330
1331     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1332         Float_t t0AC=-10000;
1333         Float_t t0A=-10000;
1334         Float_t t0C=-10000;
1335         if(flagT0T0){
1336             t0A= vevent->GetT0TOF()[1];
1337             t0C= vevent->GetT0TOF()[2];
1338             //      t0AC= vevent->GetT0TOF()[0];
1339             t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1340             resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1341             t0AC /= resT0AC*resT0AC;
1342         }
1343
1344         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1345             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1346               estimatedT0event[i]=t0AC;
1347               estimatedT0resolution[i]=resT0AC;
1348               fTOFResponse.SetT0binMask(i,6);
1349             }
1350         }
1351         else if(TMath::Abs(t0C) < t0cut){
1352             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1353               estimatedT0event[i]=t0C;
1354               estimatedT0resolution[i]=resT0C;
1355               fTOFResponse.SetT0binMask(i,4);
1356             }
1357         }
1358         else if(TMath::Abs(t0A) < t0cut){
1359             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1360               estimatedT0event[i]=t0A;
1361               estimatedT0resolution[i]=resT0A;
1362               fTOFResponse.SetT0binMask(i,2);
1363             }
1364         }
1365         else{
1366             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1367               estimatedT0event[i]=0.0;
1368               estimatedT0resolution[i]=t0spread;
1369               fTOFResponse.SetT0binMask(i,0);
1370             }
1371         }
1372         fTOFResponse.SetT0event(estimatedT0event);
1373         fTOFResponse.SetT0resolution(estimatedT0resolution);
1374     }
1375     delete [] startTime;
1376     delete [] startTimeRes;
1377     delete [] startTimeMask;
1378     delete [] estimatedT0event;
1379     delete [] estimatedT0resolution;
1380 }
1381
1382 //______________________________________________________________________________
1383 // private non cached versions of the PID calculation
1384 //
1385
1386
1387 //______________________________________________________________________________
1388 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
1389 {
1390   //
1391   // NumberOfSigmas for 'detCode'
1392   //
1393
1394   switch (detCode){
1395     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1396     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1397     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1398     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1399     default: return -999.;
1400   }
1401
1402 }
1403
1404
1405
1406 //______________________________________________________________________________
1407 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1408 {
1409   //
1410   // Calculate the number of sigmas in the ITS
1411   //
1412
1413   AliVTrack *track=(AliVTrack*)vtrack;
1414
1415   Float_t dEdx=track->GetITSsignal();
1416   if (dEdx<=0) return -999.;
1417
1418   UChar_t clumap=track->GetITSClusterMap();
1419   Int_t nPointsForPid=0;
1420   for(Int_t i=2; i<6; i++){
1421     if(clumap&(1<<i)) ++nPointsForPid;
1422   }
1423   Float_t mom=track->P();
1424
1425   //check for ITS standalone tracks
1426   Bool_t isSA=kTRUE;
1427   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1428
1429   //TODO: in case of the electron, use the SA parametrisation,
1430   //      this needs to be changed if ITS provides a parametrisation
1431   //      for electrons also for ITS+TPC tracks
1432   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1433 }
1434
1435 //______________________________________________________________________________
1436 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1437 {
1438   //
1439   // Calculate the number of sigmas in the TPC
1440   //
1441
1442   AliVTrack *track=(AliVTrack*)vtrack;
1443
1444   Double_t mom  = track->GetTPCmomentum();
1445   Double_t sig  = track->GetTPCsignal();
1446   if(fTuneMConData) sig = this->GetTPCsignalTunedOnData(track);
1447   UInt_t   sigN = track->GetTPCsignalN();
1448
1449   Double_t nSigma = -999.;
1450   if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
1451
1452   return nSigma;
1453 }
1454
1455 //______________________________________________________________________________
1456 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1457 {
1458   //
1459   // Calculate the number of sigmas in the EMCAL
1460   //
1461
1462   AliVTrack *track=(AliVTrack*)vtrack;
1463
1464   AliVCluster *matchedClus = NULL;
1465
1466   Double_t mom     = -1.;
1467   Double_t pt      = -1.;
1468   Double_t EovP    = -1.;
1469   Double_t fClsE   = -1.;
1470
1471   Int_t nMatchClus = -1;
1472   Int_t charge     = 0;
1473
1474   // Track matching
1475   nMatchClus = track->GetEMCALcluster();
1476   if(nMatchClus > -1){
1477
1478     mom    = track->P();
1479     pt     = track->Pt();
1480     charge = track->Charge();
1481
1482     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1483
1484     if(matchedClus){
1485
1486       // matched cluster is EMCAL
1487       if(matchedClus->IsEMCAL()){
1488
1489         fClsE       = matchedClus->E();
1490         EovP        = fClsE/mom;
1491
1492
1493         // NSigma value really meaningful only for electrons!
1494         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1495       }
1496     }
1497   }
1498
1499   return -999;
1500
1501 }
1502
1503
1504 //______________________________________________________________________________
1505 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1506 {
1507   //
1508   // Compute PID response of 'detCode'
1509   //
1510
1511   switch (detCode){
1512     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1513     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1514     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1515     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1516     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1517     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1518     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1519     default: return kDetNoSignal;
1520   }
1521 }
1522
1523 //______________________________________________________________________________
1524 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1525 {
1526   //
1527   // Compute PID response for the ITS
1528   //
1529
1530   // look for cached value first
1531   // only the non SA tracks are cached
1532   if (track->GetDetectorPID()){
1533     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1534   }
1535
1536   // set flat distribution (no decision)
1537   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1538
1539   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
1540     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
1541
1542   //check for ITS standalone tracks
1543   Bool_t isSA=kTRUE;
1544   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1545
1546   Double_t mom=track->P();
1547   Double_t dedx=track->GetITSsignal();
1548   Double_t momITS=mom;
1549   UChar_t clumap=track->GetITSClusterMap();
1550   Int_t nPointsForPid=0;
1551   for(Int_t i=2; i<6; i++){
1552     if(clumap&(1<<i)) ++nPointsForPid;
1553   }
1554
1555   if(nPointsForPid<3) { // track not to be used for combined PID purposes
1556     //       track->ResetStatus(AliVTrack::kITSpid);
1557     return kDetNoSignal;
1558   }
1559
1560   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1561   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
1562     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1563     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1564     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1565     //TODO: in case of the electron, use the SA parametrisation,
1566     //      this needs to be changed if ITS provides a parametrisation
1567     //      for electrons also for ITS+TPC tracks
1568     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1569     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1570       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1571     } else {
1572       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1573       mismatch=kFALSE;
1574     }
1575
1576     // Check for particles heavier than (AliPID::kSPECIES - 1)
1577     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
1578
1579   }
1580
1581   if (mismatch){
1582     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1583     return kDetNoSignal;
1584   }
1585
1586
1587   return kDetPidOk;
1588 }
1589 //______________________________________________________________________________
1590 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1591 {
1592   //
1593   // Compute PID response for the TPC
1594   //
1595
1596   // set flat distribution (no decision)
1597   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1598
1599   // check quality of the track
1600   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
1601
1602   Double_t mom = track->GetTPCmomentum();
1603
1604   Double_t dedx=track->GetTPCsignal();
1605   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1606
1607   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1608
1609   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1610     AliPID::EParticleType type=AliPID::EParticleType(j);
1611     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
1612     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
1613     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1614       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1615     } else {
1616       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1617       mismatch=kFALSE;
1618     }
1619   }
1620
1621   if (mismatch){
1622     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1623     return kDetNoSignal;
1624   }
1625
1626   return kDetPidOk;
1627 }
1628 //______________________________________________________________________________
1629 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1630 {
1631   //
1632   // Compute PID response for the
1633   //
1634
1635   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1636
1637   // set flat distribution (no decision)
1638   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1639
1640   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
1641   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
1642
1643   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
1644   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1645     AliPID::EParticleType type=AliPID::EParticleType(j);
1646     Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor;
1647
1648     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1649     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1650     if (TMath::Abs(nsigmas) > (fRange+2)) {
1651       if(nsigmas < fTOFtail)
1652         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1653       else
1654         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1655     } else{
1656       if(nsigmas < fTOFtail)
1657         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1658       else
1659         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1660     }
1661
1662     if (TMath::Abs(nsigmas)<5.){
1663       Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
1664       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
1665     }
1666   }
1667
1668   if (mismatch){
1669     return kDetMismatch;
1670   }
1671
1672   // TODO: Light nuclei
1673
1674   return kDetPidOk;
1675 }
1676 //______________________________________________________________________________
1677 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1678 {
1679   //
1680   // Compute PID response for the
1681   //
1682
1683   UInt_t TRDslicesForPID[2];
1684   SetTRDSlices(TRDslicesForPID,PIDmethod);
1685   // set flat distribution (no decision)
1686   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1687   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
1688
1689   Float_t mom[6]={0.};
1690   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
1691   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
1692   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
1693   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
1694     mom[ilayer] = track->GetTRDmomentum(ilayer);
1695     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
1696       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
1697     }
1698   }
1699   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
1700   return kDetPidOk;
1701
1702 }
1703 //______________________________________________________________________________
1704 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1705 {
1706   //
1707   // Compute PID response for the EMCAL
1708   //
1709
1710   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1711
1712   AliVCluster *matchedClus = NULL;
1713
1714   Double_t mom     = -1.;
1715   Double_t pt      = -1.;
1716   Double_t EovP    = -1.;
1717   Double_t fClsE   = -1.;
1718
1719   Int_t nMatchClus = -1;
1720   Int_t charge     = 0;
1721
1722   // Track matching
1723   nMatchClus = track->GetEMCALcluster();
1724
1725   if(nMatchClus > -1){
1726
1727     mom    = track->P();
1728     pt     = track->Pt();
1729     charge = track->Charge();
1730
1731     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1732
1733     if(matchedClus){
1734
1735       // matched cluster is EMCAL
1736       if(matchedClus->IsEMCAL()){
1737
1738         fClsE       = matchedClus->E();
1739         EovP        = fClsE/mom;
1740
1741
1742         // compute the probabilities
1743         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
1744
1745           // in case everything is OK
1746           return kDetPidOk;
1747         }
1748       }
1749     }
1750   }
1751
1752   // in all other cases set flat distribution (no decision)
1753   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
1754   return kDetNoSignal;
1755
1756 }
1757 //______________________________________________________________________________
1758 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
1759 {
1760   //
1761   // Compute PID response for the PHOS
1762   //
1763
1764   // set flat distribution (no decision)
1765   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1766   return kDetNoSignal;
1767 }
1768 //______________________________________________________________________________
1769 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1770 {
1771   //
1772   // Compute PID response for the HMPID
1773   //
1774   // set flat distribution (no decision)
1775   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1776   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
1777
1778   track->GetHMPIDpid(p);
1779
1780   return kDetPidOk;
1781 }