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 <TH2D.h>
  31 #include <TSpline.h>
  32 #include <TFile.h>
  33 #include <TArrayI.h>
  34 #include <TArrayF.h>
  35 #include <TLinearFitter.h>
  36
  37 #include <AliVEvent.h>
  38 #include <AliVTrack.h>
  39 #include <AliLog.h>
  40 #include <AliPID.h>
  41 #include <AliOADBContainer.h>
  42 #include <AliTRDPIDResponseObject.h>
  43 #include <AliTOFPIDParams.h>
  44
  45 #include "AliPIDResponse.h"
  46 #include "AliDetectorPID.h"
  47
  48 #include "AliCentrality.h"
  49
  50 ClassImp(AliPIDResponse);
  51
  52 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
  53 TNamed("PIDResponse","PIDResponse"),
  54 fITSResponse(isMC),
  55 fTPCResponse(),
  56 fTRDResponse(),
  57 fTOFResponse(),
  58 fEMCALResponse(),
  59 fRange(5.),
  60 fITSPIDmethod(kITSTruncMean),
  61 fIsMC(isMC),
  62 fCachePID(kTRUE),
  63 fOADBPath(),
  64 fCustomTPCpidResponse(),
  65 fBeamType("PP"),
  66 fLHCperiod(),
  67 fMCperiodTPC(),
  68 fMCperiodUser(),
  69 fCurrentFile(),
  70 fRecoPass(0),
  71 fRecoPassUser(-1),
  72 fRun(0),
  73 fOldRun(0),
  74 fResT0A(75.),
  75 fResT0C(65.),
  76 fResT0AC(55.),
  77 fArrPidResponseMaster(NULL),
  78 fResolutionCorrection(NULL),
  79 fOADBvoltageMaps(NULL),
  80 fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
  81 fTRDPIDResponseObject(NULL),
  82 fTOFtail(1.1),
  83 fTOFPIDParams(NULL),
  84 fEMCALPIDParams(NULL),
  85 fCurrentEvent(NULL),
  86 fCurrCentrality(0.0),
  87 fTuneMConData(kFALSE)
  88 {
  89   //
  90   // default ctor
  91   //
  92   AliLog::SetClassDebugLevel("AliPIDResponse",0);
  93   AliLog::SetClassDebugLevel("AliESDpid",0);
  94   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
  95
  96 }
  97
  98 //______________________________________________________________________________
  99 AliPIDResponse::~AliPIDResponse()
 100 {
 101   //
 102   // dtor
 103   //
 104   delete fArrPidResponseMaster;
 105   delete fTRDPIDResponseObject;
 106   delete fTOFPIDParams;
 107 }
 108
 109 //______________________________________________________________________________
 110 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
 111 TNamed(other),
 112 fITSResponse(other.fITSResponse),
 113 fTPCResponse(other.fTPCResponse),
 114 fTRDResponse(other.fTRDResponse),
 115 fTOFResponse(other.fTOFResponse),
 116 fEMCALResponse(other.fEMCALResponse),
 117 fRange(other.fRange),
 118 fITSPIDmethod(other.fITSPIDmethod),
 119 fIsMC(other.fIsMC),
 120 fCachePID(other.fCachePID),
 121 fOADBPath(other.fOADBPath),
 122 fCustomTPCpidResponse(other.fCustomTPCpidResponse),
 123 fBeamType("PP"),
 124 fLHCperiod(),
 125 fMCperiodTPC(),
 126 fMCperiodUser(other.fMCperiodUser),
 127 fCurrentFile(),
 128 fRecoPass(0),
 129 fRecoPassUser(other.fRecoPassUser),
 130 fRun(0),
 131 fOldRun(0),
 132 fResT0A(75.),
 133 fResT0C(65.),
 134 fResT0AC(55.),
 135 fArrPidResponseMaster(NULL),
 136 fResolutionCorrection(NULL),
 137 fOADBvoltageMaps(NULL),
 138 fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
 139 fTRDPIDResponseObject(NULL),
 140 fTOFtail(1.1),
 141 fTOFPIDParams(NULL),
 142 fEMCALPIDParams(NULL),
 143 fCurrentEvent(NULL),
 144 fCurrCentrality(0.0),
 145 fTuneMConData(kFALSE)
 146 {
 147   //
 148   // copy ctor
 149   //
 150 }
 151
 152 //______________________________________________________________________________
 153 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
 154 {
 155   //
 156   // copy ctor
 157   //
 158   if(this!=&other) {
 159     delete fArrPidResponseMaster;
 160     TNamed::operator=(other);
 161     fITSResponse=other.fITSResponse;
 162     fTPCResponse=other.fTPCResponse;
 163     fTRDResponse=other.fTRDResponse;
 164     fTOFResponse=other.fTOFResponse;
 165     fEMCALResponse=other.fEMCALResponse;
 166     fRange=other.fRange;
 167     fITSPIDmethod=other.fITSPIDmethod;
 168     fOADBPath=other.fOADBPath;
 169     fCustomTPCpidResponse=other.fCustomTPCpidResponse;
 170     fIsMC=other.fIsMC;
 171     fCachePID=other.fCachePID;
 172     fBeamType="PP";
 173     fLHCperiod="";
 174     fMCperiodTPC="";
 175     fMCperiodUser=other.fMCperiodUser;
 176     fCurrentFile="";
 177     fRecoPass=0;
 178     fRecoPassUser=other.fRecoPassUser;
 179     fRun=0;
 180     fOldRun=0;
 181     fResT0A=75.;
 182     fResT0C=65.;
 183     fResT0AC=55.;
 184     fArrPidResponseMaster=NULL;
 185     fResolutionCorrection=NULL;
 186     fOADBvoltageMaps=NULL;
 187         fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
 188     fTRDPIDResponseObject=NULL;
 189     fEMCALPIDParams=NULL;
 190     fTOFtail=1.1;
 191     fTOFPIDParams=NULL;
 192     fCurrentEvent=other.fCurrentEvent;
 193
 194   }
 195   return *this;
 196 }
 197
 198 //______________________________________________________________________________
 199 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
 200 {
 201   //
 202   // NumberOfSigmas for 'detCode'
 203   //
 204
 205   switch (detCode){
 206     case kDetITS: return NumberOfSigmasITS(track, type); break;
 207     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
 208     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
 209     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
 210     default: return -999.;
 211   }
 212
 213 }
 214
 215 //______________________________________________________________________________
 216 Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
 217 {
 218   //
 219   // NumberOfSigmas for 'detCode'
 220   //
 221   return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
 222 }
 223
 224 //______________________________________________________________________________
 225 // public buffered versions of the PID calculation
 226 //
 227
 228 //______________________________________________________________________________
 229 Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
 230 {
 231   //
 232   // Calculate the number of sigmas in the ITS
 233   //
 234
 235   AliVTrack *track=(AliVTrack*)vtrack;
 236
 237   // look for cached value first
 238   const AliDetectorPID *detPID=track->GetDetectorPID();
 239   const EDetector detector=kITS;
 240
 241   if ( detPID && detPID->HasNumberOfSigmas(detector)){
 242     return detPID->GetNumberOfSigmas(detector, type);
 243   } else if (fCachePID) {
 244     FillTrackDetectorPID(track, detector);
 245     detPID=track->GetDetectorPID();
 246     return detPID->GetNumberOfSigmas(detector, type);
 247   }
 248
 249   return GetNumberOfSigmasITS(track, type);
 250 }
 251
 252 //______________________________________________________________________________
 253 Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
 254 {
 255   //
 256   // Calculate the number of sigmas in the TPC
 257   //
 258
 259   AliVTrack *track=(AliVTrack*)vtrack;
 260
 261   // look for cached value first
 262   const AliDetectorPID *detPID=track->GetDetectorPID();
 263   const EDetector detector=kTPC;
 264
 265   if ( detPID && detPID->HasNumberOfSigmas(detector)){
 266     return detPID->GetNumberOfSigmas(detector, type);
 267   } else if (fCachePID) {
 268     FillTrackDetectorPID(track, detector);
 269     detPID=track->GetDetectorPID();
 270     return detPID->GetNumberOfSigmas(detector, type);
 271   }
 272
 273   return GetNumberOfSigmasTPC(track, type);
 274 }
 275
 276 //______________________________________________________________________________
 277 Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
 278                                            AliPID::EParticleType type,
 279                                            AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
 280 {
 281   //get number of sigmas according the selected TPC gain configuration scenario
 282   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
 283
 284   Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
 285
 286   return nSigma;
 287 }
 288
 289 //______________________________________________________________________________
 290 Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
 291 {
 292   //
 293   // Calculate the number of sigmas in the TOF
 294   //
 295
 296   AliVTrack *track=(AliVTrack*)vtrack;
 297
 298   // look for cached value first
 299   const AliDetectorPID *detPID=track->GetDetectorPID();
 300   const EDetector detector=kTOF;
 301
 302   if ( detPID && detPID->HasNumberOfSigmas(detector)){
 303     return detPID->GetNumberOfSigmas(detector, type);
 304   } else if (fCachePID) {
 305     FillTrackDetectorPID(track, detector);
 306     detPID=track->GetDetectorPID();
 307     return detPID->GetNumberOfSigmas(detector, type);
 308   }
 309
 310   return GetNumberOfSigmasTOF(track, type);
 311 }
 312
 313 //______________________________________________________________________________
 314 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
 315 {
 316   //
 317   // Calculate the number of sigmas in the EMCAL
 318   //
 319
 320   AliVTrack *track=(AliVTrack*)vtrack;
 321
 322   // look for cached value first
 323   const AliDetectorPID *detPID=track->GetDetectorPID();
 324   const EDetector detector=kEMCAL;
 325
 326   if ( detPID && detPID->HasNumberOfSigmas(detector)){
 327     return detPID->GetNumberOfSigmas(detector, type);
 328   } else if (fCachePID) {
 329     FillTrackDetectorPID(track, detector);
 330     detPID=track->GetDetectorPID();
 331     return detPID->GetNumberOfSigmas(detector, type);
 332   }
 333
 334   return GetNumberOfSigmasEMCAL(track, type);
 335 }
 336
 337 //______________________________________________________________________________
 338 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4])  const
 339 {
 340   //
 341   // emcal nsigma with eop and showershape
 342   //
 343   AliVTrack *track=(AliVTrack*)vtrack;
 344
 345   AliVCluster *matchedClus = NULL;
 346
 347   Double_t mom     = -1.;
 348   Double_t pt      = -1.;
 349   Double_t EovP    = -1.;
 350   Double_t fClsE   = -1.;
 351
 352   // initialize eop and shower shape parameters
 353   eop = -1.;
 354   for(Int_t i = 0; i < 4; i++){
 355     showershape[i] = -1.;
 356   }
 357
 358   Int_t nMatchClus = -1;
 359   Int_t charge     = 0;
 360
 361   // Track matching
 362   nMatchClus = track->GetEMCALcluster();
 363   if(nMatchClus > -1){
 364
 365     mom    = track->P();
 366     pt     = track->Pt();
 367     charge = track->Charge();
 368
 369     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 370
 371     if(matchedClus){
 372
 373       // matched cluster is EMCAL
 374       if(matchedClus->IsEMCAL()){
 375
 376         fClsE       = matchedClus->E();
 377         EovP        = fClsE/mom;
 378
 379         // fill used EMCAL variables here
 380         eop            = EovP; // E/p
 381         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
 382         showershape[1] = matchedClus->GetM02(); // long axis
 383         showershape[2] = matchedClus->GetM20(); // short axis
 384         showershape[3] = matchedClus->GetDispersion(); // dispersion
 385
 386         // look for cached value first
 387         const AliDetectorPID *detPID=track->GetDetectorPID();
 388         const EDetector detector=kEMCAL;
 389
 390         if ( detPID && detPID->HasNumberOfSigmas(detector)){
 391           return detPID->GetNumberOfSigmas(detector, type);
 392         } else if (fCachePID) {
 393           FillTrackDetectorPID(track, detector);
 394           detPID=track->GetDetectorPID();
 395           return detPID->GetNumberOfSigmas(detector, type);
 396         }
 397
 398         // NSigma value really meaningful only for electrons!
 399         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
 400       }
 401     }
 402   }
 403   return -999;
 404 }
 405
 406
 407 //______________________________________________________________________________
 408 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 409 {
 410   //
 411   // Compute PID response of 'detCode'
 412   //
 413
 414   switch (detCode){
 415     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
 416     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
 417     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
 418     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
 419     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
 420     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
 421     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
 422     default: return kDetNoSignal;
 423   }
 424 }
 425
 426 //______________________________________________________________________________
 427 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 428 {
 429   //
 430   // Compute PID response of 'detCode'
 431   //
 432
 433   return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
 434 }
 435
 436 //______________________________________________________________________________
 437 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 438 {
 439   //
 440   // Compute PID response for the ITS
 441   //
 442
 443   // look for cached value first
 444   const AliDetectorPID *detPID=track->GetDetectorPID();
 445   const EDetector detector=kITS;
 446
 447   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 448     return detPID->GetRawProbability(detector, p, nSpecies);
 449   } else if (fCachePID) {
 450     FillTrackDetectorPID(track, detector);
 451     detPID=track->GetDetectorPID();
 452     return detPID->GetRawProbability(detector, p, nSpecies);
 453   }
 454
 455   return GetComputeITSProbability(track, nSpecies, p);
 456 }
 457 //______________________________________________________________________________
 458 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 459 {
 460   //
 461   // Compute PID response for the TPC
 462   //
 463
 464   // look for cached value first
 465   const AliDetectorPID *detPID=track->GetDetectorPID();
 466   const EDetector detector=kTPC;
 467
 468   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 469     return detPID->GetRawProbability(detector, p, nSpecies);
 470   } else if (fCachePID) {
 471     FillTrackDetectorPID(track, detector);
 472     detPID=track->GetDetectorPID();
 473     return detPID->GetRawProbability(detector, p, nSpecies);
 474   }
 475
 476   return GetComputeTPCProbability(track, nSpecies, p);
 477 }
 478 //______________________________________________________________________________
 479 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 480 {
 481   //
 482   // Compute PID response for the
 483   //
 484
 485   const AliDetectorPID *detPID=track->GetDetectorPID();
 486   const EDetector detector=kTOF;
 487
 488   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 489     return detPID->GetRawProbability(detector, p, nSpecies);
 490   } else if (fCachePID) {
 491     FillTrackDetectorPID(track, detector);
 492     detPID=track->GetDetectorPID();
 493     return detPID->GetRawProbability(detector, p, nSpecies);
 494   }
 495
 496   return GetComputeTOFProbability(track, nSpecies, p);
 497 }
 498 //______________________________________________________________________________
 499 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
 500 {
 501   //
 502   // Compute PID response for the
 503   //
 504
 505   const AliDetectorPID *detPID=track->GetDetectorPID();
 506   const EDetector detector=kTRD;
 507
 508   // chacke only for the default method (1d at the moment)
 509   if (PIDmethod!=AliTRDPIDResponse::kLQ1D) return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
 510
 511   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 512     return detPID->GetRawProbability(detector, p, nSpecies);
 513   } else if (fCachePID) {
 514     FillTrackDetectorPID(track, detector);
 515     detPID=track->GetDetectorPID();
 516     return detPID->GetRawProbability(detector, p, nSpecies);
 517   }
 518
 519   return GetComputeTRDProbability(track, nSpecies, p);
 520 }
 521 //______________________________________________________________________________
 522 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 523 {
 524   //
 525   // Compute PID response for the EMCAL
 526   //
 527
 528   const AliDetectorPID *detPID=track->GetDetectorPID();
 529   const EDetector detector=kEMCAL;
 530
 531   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 532     return detPID->GetRawProbability(detector, p, nSpecies);
 533   } else if (fCachePID) {
 534     FillTrackDetectorPID(track, detector);
 535     detPID=track->GetDetectorPID();
 536     return detPID->GetRawProbability(detector, p, nSpecies);
 537   }
 538
 539   return GetComputeEMCALProbability(track, nSpecies, p);
 540 }
 541 //______________________________________________________________________________
 542 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
 543 {
 544   //
 545   // Compute PID response for the PHOS
 546   //
 547
 548   // look for cached value first
 549 //   if (track->GetDetectorPID()){
 550 //     return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
 551 //   }
 552
 553   // set flat distribution (no decision)
 554   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 555   return kDetNoSignal;
 556 }
 557 //______________________________________________________________________________
 558 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 559 {
 560   //
 561   // Compute PID response for the HMPID
 562   //
 563
 564   const AliDetectorPID *detPID=track->GetDetectorPID();
 565   const EDetector detector=kHMPID;
 566
 567   if ( detPID && detPID->HasRawProbabilitiy(detector)){
 568     return detPID->GetRawProbability(detector, p, nSpecies);
 569   } else if (fCachePID) {
 570     FillTrackDetectorPID(track, detector);
 571     detPID=track->GetDetectorPID();
 572     return detPID->GetRawProbability(detector, p, nSpecies);
 573   }
 574
 575   return GetComputeHMPIDProbability(track, nSpecies, p);
 576 }
 577
 578 //______________________________________________________________________________
 579 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
 580 {
 581   //
 582   // Apply settings for the current event
 583   //
 584   fRecoPass=pass;
 585
 586
 587   fCurrentEvent=NULL;
 588   if (!event) return;
 589   fCurrentEvent=event;
 590   if (run>0) fRun=run;
 591   else fRun=event->GetRunNumber();
 592
 593   if (fRun!=fOldRun){
 594     ExecNewRun();
 595     fOldRun=fRun;
 596   }
 597
 598   //TPC resolution parametrisation PbPb
 599   if ( fResolutionCorrection ){
 600     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
 601     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
 602   }
 603
 604   //TOF resolution
 605   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 606
 607
 608   // Get and set centrality
 609   AliCentrality *centrality = event->GetCentrality();
 610   if(centrality){
 611     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
 612   }
 613   else{
 614     fCurrCentrality = -1;
 615   }
 616 }
 617
 618 //______________________________________________________________________________
 619 void AliPIDResponse::ExecNewRun()
 620 {
 621   //
 622   // Things to Execute upon a new run
 623   //
 624   SetRecoInfo();
 625
 626   SetITSParametrisation();
 627
 628   SetTPCPidResponseMaster();
 629   SetTPCParametrisation();
 630   SetTPCEtaMaps();
 631
 632   SetTRDPidResponseMaster();
 633   InitializeTRDResponse();
 634
 635   SetEMCALPidResponseMaster();
 636   InitializeEMCALResponse();
 637
 638   SetTOFPidResponseMaster();
 639   InitializeTOFResponse();
 640
 641   if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
 642 }
 643
 644 //______________________________________________________________________________
 645 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
 646 {
 647   //
 648   // Get TPC multiplicity in bins of 150
 649   //
 650
 651   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
 652   Double_t tpcMulti=0.;
 653   if(vertexTPC){
 654     Double_t vertexContribTPC=vertexTPC->GetNContributors();
 655     tpcMulti=vertexContribTPC/150.;
 656     if (tpcMulti>20.) tpcMulti=20.;
 657   }
 658
 659   return tpcMulti;
 660 }
 661
 662 //______________________________________________________________________________
 663 void AliPIDResponse::SetRecoInfo()
 664 {
 665   //
 666   // Set reconstruction information
 667   //
 668
 669   //reset information
 670   fLHCperiod="";
 671   fMCperiodTPC="";
 672
 673   fBeamType="";
 674
 675   fBeamType="PP";
 676
 677   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
 678   TPRegexp reg12a17("LHC1[2-3][a-z]");
 679
 680   //find the period by run number (UGLY, but not stored in ESD and AOD... )
 681   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
 682   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
 683   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
 684   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 685   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 686   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 687   else if (fRun>=136851&&fRun<=139846) {
 688     fLHCperiod="LHC10H";
 689     fMCperiodTPC="LHC10H8";
 690     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 691     fBeamType="PBPB";
 692   }
 693   else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 694   //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
 695   else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 696   else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 697   // also for 11e (159650-162750),f(162751-165771) use 11d
 698   else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 699   else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 700
 701   else if (fRun>=165772 && fRun<=170718) {
 702     fLHCperiod="LHC11H";
 703     fMCperiodTPC="LHC11A10";
 704     fBeamType="PBPB";
 705     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
 706   }
 707   if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 708   // for the moment use LHC12b parameters up to LHC12e
 709   if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 710 //   if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 711 //   if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 712 //   if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 713
 714 //   if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 715 //   if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 716 //   if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
 717 // for the moment use 12g parametrisation for all full gain runs (LHC12f+)
 718   if (fRun >= 186636  ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
 719
 720   //exception new pp MC productions from 2011
 721   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
 722   // exception for 11f1
 723   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 724   // exception for 12f1a, 12f1b and 12i3
 725   if (fCurrentFile.Contains("LHC12f1a/") || fCurrentFile.Contains("LHC12f1b/")
 726       || fCurrentFile.Contains("LHC12i3/")) fMCperiodTPC="LHC12F1";
 727 }
 728
 729 //______________________________________________________________________________
 730 void AliPIDResponse::SetITSParametrisation()
 731 {
 732   //
 733   // Set the ITS parametrisation
 734   //
 735 }
 736
 737
 738 //______________________________________________________________________________
 739 void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
 740 {
 741   if (h->GetBinContent(binX, binY) <= 1e-4)
 742     return; // Reject bins without content (within some numerical precision) or with strange content
 743
 744   Double_t coord[2] = {0, 0};
 745   coord[0] = h->GetXaxis()->GetBinCenter(binX);
 746   coord[1] = h->GetYaxis()->GetBinCenter(binY);
 747   Double_t binError = h->GetBinError(binX, binY);
 748   if (binError <= 0) {
 749     binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
 750     printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
 751   }
 752   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
 753 }
 754
 755
 756 //______________________________________________________________________________
 757 TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
 758 {
 759   if (!h)
 760     return 0x0;
 761
 762   // Interpolate to finer map
 763   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
 764
 765   Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
 766   Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
 767   Int_t nBinsX = 30;
 768   // Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
 769   // scale the number of bins correspondingly
 770   Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
 771   Int_t nBinsXrefined = nBinsX * refineFactorX;
 772   Int_t nBinsYrefined = nBinsY * refineFactorY;
 773
 774   TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()),  Form("%s (refined)", h->GetTitle()),
 775                             nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
 776                             nBinsYrefined, lowerMapBoundY, upperMapBoundY);
 777
 778   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
 779     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
 780
 781       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
 782
 783       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
 784       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
 785
 786       /*OLD
 787       linExtrapolation->ClearPoints();
 788
 789       // For interpolation: Just take the corresponding bin from the old histo.
 790       // For extrapolation: take the last available bin from the old histo.
 791       // If the boundaries are to be skipped, also skip the corresponding bins
 792       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
 793       if (oldBinX < 1)
 794         oldBinX = 1;
 795       if (oldBinX > nBinsX)
 796         oldBinX = nBinsX;
 797
 798       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
 799       if (oldBinY < 1)
 800         oldBinY = 1;
 801       if (oldBinY > nBinsY)
 802         oldBinY = nBinsY;
 803
 804       // Neighbours left column
 805       if (oldBinX >= 2) {
 806         if (oldBinY >= 2) {
 807           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
 808         }
 809
 810         AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
 811
 812         if (oldBinY < nBinsY) {
 813           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
 814         }
 815       }
 816
 817       // Neighbours (and point itself) same column
 818       if (oldBinY >= 2) {
 819         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
 820       }
 821
 822       AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
 823
 824       if (oldBinY < nBinsY) {
 825         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
 826       }
 827
 828       // Neighbours right column
 829       if (oldBinX < nBinsX) {
 830         if (oldBinY >= 2) {
 831           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
 832         }
 833
 834         AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
 835
 836         if (oldBinY < nBinsY) {
 837           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
 838         }
 839       }
 840
 841
 842       // Fit 2D-hyperplane
 843       if (linExtrapolation->GetNpoints() <= 0)
 844         continue;
 845
 846       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
 847         continue;
 848
 849       // Fill the bin of the refined histogram with the extrapolated value
 850       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
 851                                  + linExtrapolation->GetParameter(2) * centerY;
 852       */
 853       Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
 854       hRefined->SetBinContent(binX, binY, interpolatedValue);
 855     }
 856   }
 857
 858
 859   // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
 860   // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
 861   // Assume line through these points and extropolate to last bin of refined map
 862   const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
 863   const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
 864
 865   const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
 866
 867   const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
 868   const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
 869
 870   for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
 871     Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
 872
 873     const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
 874     const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
 875
 876     const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
 877     const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
 878
 879
 880     const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
 881     const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
 882
 883     const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
 884     const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
 885
 886     for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
 887       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
 888
 889       if (centerX < firstOldXbinCenter) {
 890         Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
 891         hRefined->SetBinContent(binX, binY, extrapolatedValue);
 892       }
 893       else if (centerX <= lastOldXbinCenter) {
 894         continue;
 895       }
 896       else {
 897         Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
 898         hRefined->SetBinContent(binX, binY, extrapolatedValue);
 899       }
 900     }
 901   }
 902
 903   delete linExtrapolation;
 904
 905   return hRefined;
 906 }
 907
 908 //______________________________________________________________________________
 909 void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
 910                                    Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
 911 {
 912   //
 913   // Load the TPC eta correction maps from the OADB
 914   //
 915
 916   if (fUseTPCEtaCorrection == kFALSE) {
 917     // Disable eta correction via setting no maps
 918     if (!fTPCResponse.SetEtaCorrMap(0x0))
 919       AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
 920     else
 921       AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
 922
 923     if (!fTPCResponse.SetSigmaParams(0x0, 0))
 924       AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
 925     else
 926       AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
 927
 928     return;
 929   }
 930
 931   TString dataType = "DATA";
 932   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
 933
 934   if (fIsMC)  {
 935     if (!fTuneMConData) {
 936       period=fMCperiodTPC;
 937       dataType="MC";
 938     }
 939     fRecoPass = 1;
 940
 941     if (!fTuneMConData && fMCperiodTPC.IsNull()) {
 942       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
 943       return;
 944     }
 945   }
 946
 947   Int_t recopass = fRecoPass;
 948   if (fTuneMConData)
 949     recopass = fRecoPassUser;
 950
 951   TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
 952
 953   AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
 954
 955   // Invalidate old maps
 956   fTPCResponse.SetEtaCorrMap(0x0);
 957   fTPCResponse.SetSigmaParams(0x0, 0);
 958
 959   // Load the eta correction maps
 960   AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
 961
 962   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
 963                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
 964   if (statusCont) {
 965     AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
 966   }
 967   else {
 968     AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
 969
 970     TH2D* etaMap = 0x0;
 971
 972     if (fIsMC && !fTuneMConData) {
 973       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
 974       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
 975       if (!etaMap) {
 976         // Try default object
 977         etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
 978       }
 979     }
 980     else {
 981       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
 982     }
 983
 984
 985     if (!etaMap) {
 986       AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
 987     }
 988     else {
 989       TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
 990
 991       if (etaMapRefined) {
 992         if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
 993           AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
 994           fTPCResponse.SetEtaCorrMap(0x0);
 995         }
 996         else {
 997           AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
 998                        refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle()));
 999         }
1000
1001         delete etaMapRefined;
1002       }
1003       else {
1004         AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
1005       }
1006     }
1007   }
1008
1009   // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
1010   AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
1011
1012   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
1013                                              Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
1014   if (statusCont) {
1015     AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
1016   }
1017   else {
1018     AliInfo(Form("Loading TPC eta sigma map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
1019
1020     TObjArray* etaSigmaPars = 0x0;
1021
1022     if (fIsMC && !fTuneMConData) {
1023       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
1024       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
1025       if (!etaSigmaPars) {
1026         // Try default object
1027         etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
1028       }
1029     }
1030     else {
1031       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
1032     }
1033
1034     if (!etaSigmaPars) {
1035       AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
1036     }
1037     else {
1038       TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
1039       TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
1040       Double_t sigmaPar0 = 0.0;
1041
1042       if (sigmaPar0Info) {
1043         TString sigmaPar0String = sigmaPar0Info->GetTitle();
1044         sigmaPar0 = sigmaPar0String.Atof();
1045       }
1046       else {
1047         // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
1048         etaSigmaPar1Map = 0x0;
1049       }
1050
1051       TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
1052
1053
1054       if (etaSigmaPar1MapRefined) {
1055         if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
1056           AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
1057           fTPCResponse.SetSigmaParams(0x0, 0);
1058         }
1059         else {
1060           AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
1061                        refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle()));
1062         }
1063
1064         delete etaSigmaPar1MapRefined;
1065       }
1066       else {
1067         AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
1068                       fRun));
1069       }
1070     }
1071   }
1072 }
1073
1074 //______________________________________________________________________________
1075 void AliPIDResponse::SetTPCPidResponseMaster()
1076 {
1077   //
1078   // Load the TPC pid response functions from the OADB
1079   // Load the TPC voltage maps from OADB
1080   //
1081   //don't load twice for the moment
1082    if (fArrPidResponseMaster) return;
1083
1084
1085   //reset the PID response functions
1086   delete fArrPidResponseMaster;
1087   fArrPidResponseMaster=NULL;
1088
1089   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
1090   TFile *f=NULL;
1091   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
1092
1093   TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
1094   f=TFile::Open(fileNamePIDresponse.Data());
1095   if (f && f->IsOpen() && !f->IsZombie()){
1096     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
1097   }
1098   delete f;
1099
1100   TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
1101   f=TFile::Open(fileNameVoltageMaps.Data());
1102   if (f && f->IsOpen() && !f->IsZombie()){
1103     fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
1104   }
1105   delete f;
1106
1107   if (!fArrPidResponseMaster){
1108     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
1109     return;
1110   }
1111   fArrPidResponseMaster->SetOwner();
1112
1113   if (!fOADBvoltageMaps)
1114   {
1115     AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
1116   }
1117   fArrPidResponseMaster->SetOwner();
1118 }
1119
1120 //______________________________________________________________________________
1121 void AliPIDResponse::SetTPCParametrisation()
1122 {
1123   //
1124   // Change BB parametrisation for current run
1125   //
1126
1127   //
1128   //reset old splines
1129   //
1130   fTPCResponse.ResetSplines();
1131
1132   if (fLHCperiod.IsNull()) {
1133     AliError("No period set, not changing parametrisation");
1134     return;
1135   }
1136
1137   //
1138   // Set default parametrisations for data and MC
1139   //
1140
1141   //data type
1142   TString datatype="DATA";
1143   //in case of mc fRecoPass is per default 1
1144   if (fIsMC) {
1145       if(!fTuneMConData) datatype="MC";
1146       fRecoPass=1;
1147   }
1148
1149   // period
1150   TString period=fLHCperiod;
1151   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
1152
1153   Int_t recopass = fRecoPass;
1154   if(fTuneMConData) recopass = fRecoPassUser;
1155
1156   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1157   Bool_t found=kFALSE;
1158   //
1159   //set the new PID splines
1160   //
1161   if (fArrPidResponseMaster){
1162     //for MC don't use period information
1163     //if (fIsMC) period="[A-Z0-9]*";
1164     //for MC use MC period information
1165     //pattern for the default entry (valid for all particles)
1166     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1167
1168     //find particle id and gain scenario
1169     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1170     {
1171       TObject *grAll=NULL;
1172       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1173       gainScenario.ToUpper();
1174       //loop over entries and filter them
1175       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1176       {
1177         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1178         if (responseFunction==NULL) continue;
1179         TString responseName=responseFunction->GetName();
1180
1181         if (!reg.MatchB(responseName)) continue;
1182
1183         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1184         TObject* tmp=NULL;
1185         tmp=arr->At(1); if (!tmp) continue;
1186         TString particleName=tmp->GetName();
1187         tmp=arr->At(3); if (!tmp) continue;
1188         TString gainScenarioName=tmp->GetName();
1189         delete arr;
1190         if (particleName.IsNull()) continue;
1191         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1192         else
1193         {
1194           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1195           {
1196             TString particle=AliPID::ParticleName(ispec);
1197             particle.ToUpper();
1198             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1199             if ( particle == particleName && gainScenario == gainScenarioName )
1200             {
1201               fTPCResponse.SetResponseFunction( responseFunction,
1202                                                 (AliPID::EParticleType)ispec,
1203                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1204               fTPCResponse.SetUseDatabase(kTRUE);
1205               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
1206               found=kTRUE;
1207               break;
1208             }
1209           }
1210         }
1211       }
1212
1213       // Retrieve responsefunction for pions - will (if available) be used for muons if there are no dedicated muon splines.
1214       // For light nuclei, try to set the proton spline, if no dedicated splines are available.
1215       // In both cases: Use default splines, if no dedicated splines and no pion/proton splines are available.
1216       TObject* responseFunctionPion = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kPion,
1217                                                                         (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1218       TObject* responseFunctionProton = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kProton,
1219                                                                           (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1220
1221       for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1222       {
1223         if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1224           (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1225         {
1226           if (ispec == AliPID::kMuon) { // Muons
1227             if (responseFunctionPion) {
1228               fTPCResponse.SetResponseFunction( responseFunctionPion,
1229                                                 (AliPID::EParticleType)ispec,
1230                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1231               fTPCResponse.SetUseDatabase(kTRUE);
1232               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionPion->GetName()));
1233               found=kTRUE;
1234             }
1235             else if (grAll) {
1236               fTPCResponse.SetResponseFunction( grAll,
1237                                                 (AliPID::EParticleType)ispec,
1238                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1239               fTPCResponse.SetUseDatabase(kTRUE);
1240               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1241               found=kTRUE;
1242             }
1243             //else
1244             //  AliError(Form("No splines found for muons (also no pion splines and no default splines) for gain scenario %d!", igainScenario));
1245           }
1246           else if (ispec >= AliPID::kSPECIES) { // Light nuclei
1247             if (responseFunctionProton) {
1248               fTPCResponse.SetResponseFunction( responseFunctionProton,
1249                                                 (AliPID::EParticleType)ispec,
1250                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1251               fTPCResponse.SetUseDatabase(kTRUE);
1252               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionProton->GetName()));
1253               found=kTRUE;
1254             }
1255             else if (grAll) {
1256               fTPCResponse.SetResponseFunction( grAll,
1257                                                 (AliPID::EParticleType)ispec,
1258                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1259               fTPCResponse.SetUseDatabase(kTRUE);
1260               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1261               found=kTRUE;
1262             }
1263             //else
1264             //  AliError(Form("No splines found for species %d (also no proton splines and no default splines) for gain scenario %d!",
1265             //                ispec, igainScenario));
1266           }
1267         }
1268       }
1269     }
1270   }
1271   else AliInfo("no fArrPidResponseMaster");
1272
1273   if (!found){
1274     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1275   }
1276
1277   //
1278   // Setup resolution parametrisation
1279   //
1280
1281   //default
1282   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1283
1284   if (fRun>=122195){
1285     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1286   }
1287
1288   if (fRun>=186636){
1289 //   if (fRun>=188356){
1290     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1291   }
1292
1293   if (fArrPidResponseMaster)
1294   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1295
1296   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1297
1298   //read in the voltage map
1299   TVectorF* gsm = 0x0;
1300   if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1301   if (gsm)
1302   {
1303     fTPCResponse.SetVoltageMap(*gsm);
1304     TString vals;
1305     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1306     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1307     AliInfo(vals.Data());
1308     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1309     AliInfo(vals.Data());
1310     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1311     AliInfo(vals.Data());
1312     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1313     AliInfo(vals.Data());
1314   }
1315   else AliInfo("no voltage map, ideal default assumed");
1316 }
1317
1318 //______________________________________________________________________________
1319 void AliPIDResponse::SetTRDPidResponseMaster()
1320 {
1321   //
1322   // Load the TRD pid params and references from the OADB
1323   //
1324   if(fTRDPIDResponseObject) return;
1325   AliOADBContainer contParams("contParams");
1326
1327   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1328   if(statusResponse){
1329     AliError("Failed initializing PID Response Object from OADB");
1330   } else {
1331     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1332     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1333     if(!fTRDPIDResponseObject){
1334       AliError(Form("TRD Response not found in run %d", fRun));
1335     }
1336   }
1337 }
1338
1339 //______________________________________________________________________________
1340 void AliPIDResponse::InitializeTRDResponse(){
1341   //
1342   // Set PID Params and references to the TRD PID response
1343   //
1344     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1345 }
1346
1347 //______________________________________________________________________________
1348 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1349
1350     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1351         // backward compatibility for setting with 8 slices
1352         TRDslicesForPID[0] = 0;
1353         TRDslicesForPID[1] = 7;
1354     }
1355     else{
1356         if(method==AliTRDPIDResponse::kLQ1D){
1357             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1358             TRDslicesForPID[1] = 0;
1359         }
1360         if(method==AliTRDPIDResponse::kLQ2D){
1361             TRDslicesForPID[0] = 1;
1362             TRDslicesForPID[1] = 7;
1363         }
1364     }
1365     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1366 }
1367
1368 //______________________________________________________________________________
1369 void AliPIDResponse::SetTOFPidResponseMaster()
1370 {
1371   //
1372   // Load the TOF pid params from the OADB
1373   //
1374
1375   if (fTOFPIDParams) delete fTOFPIDParams;
1376   fTOFPIDParams=NULL;
1377
1378   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1379   if (oadbf && oadbf->IsOpen()) {
1380     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1381     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1382     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1383     oadbf->Close();
1384     delete oadbc;
1385   }
1386   delete oadbf;
1387
1388   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1389 }
1390
1391 //______________________________________________________________________________
1392 void AliPIDResponse::InitializeTOFResponse(){
1393   //
1394   // Set PID Params to the TOF PID response
1395   //
1396
1397   AliInfo("TOF PID Params loaded from OADB");
1398   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1399   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1400   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1401                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1402
1403   for (Int_t i=0;i<4;i++) {
1404     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1405   }
1406   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1407
1408   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1409   Float_t t0Spread[4];
1410   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1411   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]));
1412   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1413   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1414   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1415     fResT0AC=t0Spread[3];
1416     fResT0A=TMath::Sqrt(a);
1417     fResT0C=TMath::Sqrt(c);
1418   } else {
1419     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1420     fResT0A=75.;
1421     fResT0C=65.;
1422     fResT0AC=55.;
1423   }
1424   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1425
1426 }
1427
1428
1429 //______________________________________________________________________________
1430 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1431   //
1432   // Check whether track is identified as electron under a given electron efficiency hypothesis
1433     //
1434
1435   Double_t probs[AliPID::kSPECIES];
1436   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1437
1438   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1439   // Take mean of the TRD momenta in the given tracklets
1440   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1441   Int_t nmomenta = 0;
1442   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1443     if(vtrack->GetTRDmomentum(iPl) > 0.){
1444       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1445     }
1446   }
1447   p = TMath::Mean(nmomenta, trdmomenta);
1448
1449   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1450 }
1451
1452 //______________________________________________________________________________
1453 void AliPIDResponse::SetEMCALPidResponseMaster()
1454 {
1455   //
1456   // Load the EMCAL pid response functions from the OADB
1457   //
1458   TObjArray* fEMCALPIDParamsRun      = NULL;
1459   TObjArray* fEMCALPIDParamsPass     = NULL;
1460
1461   if(fEMCALPIDParams) return;
1462   AliOADBContainer contParams("contParams");
1463
1464   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1465   if(statusPars){
1466     AliError("Failed initializing PID Params from OADB");
1467   }
1468   else {
1469     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1470
1471     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1472     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1473     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1474
1475     if(!fEMCALPIDParams){
1476       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1477       AliInfo("Will take the standard LHC11d instead ...");
1478
1479       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1480       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1481       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1482
1483       if(!fEMCALPIDParams){
1484         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1485       }
1486     }
1487   }
1488 }
1489
1490 //______________________________________________________________________________
1491 void AliPIDResponse::InitializeEMCALResponse(){
1492   //
1493   // Set PID Params to the EMCAL PID response
1494   //
1495   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1496
1497 }
1498
1499 //______________________________________________________________________________
1500 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1501 {
1502   //
1503   // create detector PID information and setup the transient pointer in the track
1504   //
1505
1506   // check if detector number is inside accepted range
1507   if (detector == kNdetectors) return;
1508
1509   // get detector pid
1510   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1511   if (!detPID) {
1512     detPID=new AliDetectorPID;
1513     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1514   }
1515
1516   //check if values exist
1517   if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return;
1518
1519   //TODO: which particles to include? See also the loops below...
1520   Double_t values[AliPID::kSPECIESC]={0};
1521
1522   //nsigmas
1523   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1524     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1525   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC);
1526
1527   //probabilities
1528   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1529   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1530 }
1531
1532 //______________________________________________________________________________
1533 void AliPIDResponse::FillTrackDetectorPID()
1534 {
1535   //
1536   // create detector PID information and setup the transient pointer in the track
1537   //
1538
1539   if (!fCurrentEvent) return;
1540
1541   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1542     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1543     if (!track) continue;
1544
1545     for (Int_t idet=0; idet<kNdetectors; ++idet){
1546       FillTrackDetectorPID(track, (EDetector)idet);
1547     }
1548   }
1549 }
1550
1551 //______________________________________________________________________________
1552 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1553   //
1554   // Set TOF response function
1555   // Input option for event_time used
1556   //
1557
1558     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1559     if(t0spread < 10) t0spread = 80;
1560
1561     // T0 from TOF algorithm
1562
1563     Bool_t flagT0TOF=kFALSE;
1564     Bool_t flagT0T0=kFALSE;
1565     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1566     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1567     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1568
1569     // T0-TOF arrays
1570     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1571     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1572     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1573       estimatedT0event[i]=0.0;
1574       estimatedT0resolution[i]=0.0;
1575       startTimeMask[i] = 0;
1576     }
1577
1578     Float_t resT0A=fResT0A;
1579     Float_t resT0C=fResT0C;
1580     Float_t resT0AC=fResT0AC;
1581     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1582         flagT0T0=kTRUE;
1583     }
1584
1585
1586     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1587
1588     if (tofHeader) { // read global info and T0-TOF
1589       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1590       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1591       if(t0spread < 10) t0spread = 80;
1592
1593       flagT0TOF=kTRUE;
1594       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1595         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1596         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1597         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1598       }
1599
1600       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1601       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1602       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1603       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1604         Int_t icurrent = (Int_t)ibin->GetAt(j);
1605         startTime[icurrent]=t0Bin->GetAt(j);
1606         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1607         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1608       }
1609     }
1610
1611     // for cut of 3 sigma on t0 spread
1612     Float_t t0cut = 3 * t0spread;
1613     if(t0cut < 500) t0cut = 500;
1614
1615     if(option == kFILL_T0){ // T0-FILL is used
1616         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1617           estimatedT0event[i]=0.0;
1618           estimatedT0resolution[i]=t0spread;
1619         }
1620         fTOFResponse.SetT0event(estimatedT0event);
1621         fTOFResponse.SetT0resolution(estimatedT0resolution);
1622     }
1623
1624     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1625         if(flagT0TOF){
1626             fTOFResponse.SetT0event(startTime);
1627             fTOFResponse.SetT0resolution(startTimeRes);
1628             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1629               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1630               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1631             }
1632         }
1633         else{
1634             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1635               estimatedT0event[i]=0.0;
1636               estimatedT0resolution[i]=t0spread;
1637               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1638             }
1639             fTOFResponse.SetT0event(estimatedT0event);
1640             fTOFResponse.SetT0resolution(estimatedT0resolution);
1641         }
1642     }
1643     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1644         Float_t t0AC=-10000;
1645         Float_t t0A=-10000;
1646         Float_t t0C=-10000;
1647         if(flagT0T0){
1648             t0A= vevent->GetT0TOF()[1];
1649             t0C= vevent->GetT0TOF()[2];
1650         //      t0AC= vevent->GetT0TOF()[0];
1651         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1652         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1653         t0AC /= resT0AC*resT0AC;
1654         }
1655
1656         Float_t t0t0Best = 0;
1657         Float_t t0t0BestRes = 9999;
1658         Int_t t0used=0;
1659         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1660             t0t0Best = t0AC;
1661             t0t0BestRes = resT0AC;
1662             t0used=6;
1663         }
1664         else if(TMath::Abs(t0C) < t0cut){
1665             t0t0Best = t0C;
1666             t0t0BestRes = resT0C;
1667             t0used=4;
1668         }
1669         else if(TMath::Abs(t0A) < t0cut){
1670             t0t0Best = t0A;
1671             t0t0BestRes = resT0A;
1672             t0used=2;
1673         }
1674
1675         if(flagT0TOF){ // if T0-TOF info is available
1676             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1677                 if(t0t0BestRes < 999){
1678                   if(startTimeRes[i] < t0spread){
1679                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1680                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1681                     estimatedT0event[i]=t0best / wtot;
1682                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1683                     startTimeMask[i] = t0used+1;
1684                   }
1685                   else {
1686                     estimatedT0event[i]=t0t0Best;
1687                     estimatedT0resolution[i]=t0t0BestRes;
1688                     startTimeMask[i] = t0used;
1689                   }
1690                 }
1691                 else{
1692                   estimatedT0event[i]=startTime[i];
1693                   estimatedT0resolution[i]=startTimeRes[i];
1694                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1695                 }
1696                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1697             }
1698             fTOFResponse.SetT0event(estimatedT0event);
1699             fTOFResponse.SetT0resolution(estimatedT0resolution);
1700         }
1701         else{ // if no T0-TOF info is available
1702             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1703               fTOFResponse.SetT0binMask(i,t0used);
1704               if(t0t0BestRes < 999){
1705                 estimatedT0event[i]=t0t0Best;
1706                 estimatedT0resolution[i]=t0t0BestRes;
1707               }
1708               else{
1709                 estimatedT0event[i]=0.0;
1710                 estimatedT0resolution[i]=t0spread;
1711               }
1712             }
1713             fTOFResponse.SetT0event(estimatedT0event);
1714             fTOFResponse.SetT0resolution(estimatedT0resolution);
1715         }
1716     }
1717
1718     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1719         Float_t t0AC=-10000;
1720         Float_t t0A=-10000;
1721         Float_t t0C=-10000;
1722         if(flagT0T0){
1723             t0A= vevent->GetT0TOF()[1];
1724             t0C= vevent->GetT0TOF()[2];
1725         //      t0AC= vevent->GetT0TOF()[0];
1726         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1727         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1728         t0AC /= resT0AC*resT0AC;
1729         }
1730
1731         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1732             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1733               estimatedT0event[i]=t0AC;
1734               estimatedT0resolution[i]=resT0AC;
1735               fTOFResponse.SetT0binMask(i,6);
1736             }
1737         }
1738         else if(TMath::Abs(t0C) < t0cut){
1739             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1740               estimatedT0event[i]=t0C;
1741               estimatedT0resolution[i]=resT0C;
1742               fTOFResponse.SetT0binMask(i,4);
1743             }
1744         }
1745         else if(TMath::Abs(t0A) < t0cut){
1746             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1747               estimatedT0event[i]=t0A;
1748               estimatedT0resolution[i]=resT0A;
1749               fTOFResponse.SetT0binMask(i,2);
1750             }
1751         }
1752         else{
1753             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1754               estimatedT0event[i]=0.0;
1755               estimatedT0resolution[i]=t0spread;
1756               fTOFResponse.SetT0binMask(i,0);
1757             }
1758         }
1759         fTOFResponse.SetT0event(estimatedT0event);
1760         fTOFResponse.SetT0resolution(estimatedT0resolution);
1761     }
1762     delete [] startTime;
1763     delete [] startTimeRes;
1764     delete [] startTimeMask;
1765     delete [] estimatedT0event;
1766     delete [] estimatedT0resolution;
1767 }
1768
1769 //______________________________________________________________________________
1770 // private non cached versions of the PID calculation
1771 //
1772
1773
1774 //______________________________________________________________________________
1775 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
1776 {
1777   //
1778   // NumberOfSigmas for 'detCode'
1779   //
1780
1781   switch (detCode){
1782     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1783     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1784     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1785     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1786     default: return -999.;
1787   }
1788
1789 }
1790
1791
1792
1793 //______________________________________________________________________________
1794 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1795 {
1796   //
1797   // Calculate the number of sigmas in the ITS
1798   //
1799
1800   AliVTrack *track=(AliVTrack*)vtrack;
1801
1802   Float_t dEdx=track->GetITSsignal();
1803   if (dEdx<=0) return -999.;
1804
1805   UChar_t clumap=track->GetITSClusterMap();
1806   Int_t nPointsForPid=0;
1807   for(Int_t i=2; i<6; i++){
1808     if(clumap&(1<<i)) ++nPointsForPid;
1809   }
1810   Float_t mom=track->P();
1811
1812   //check for ITS standalone tracks
1813   Bool_t isSA=kTRUE;
1814   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1815
1816   //TODO: in case of the electron, use the SA parametrisation,
1817   //      this needs to be changed if ITS provides a parametrisation
1818   //      for electrons also for ITS+TPC tracks
1819   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1820 }
1821
1822 //______________________________________________________________________________
1823 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1824 {
1825   //
1826   // Calculate the number of sigmas in the TPC
1827   //
1828
1829   AliVTrack *track=(AliVTrack*)vtrack;
1830
1831   Double_t nSigma = -999.;
1832
1833   if (fTuneMConData)
1834     this->GetTPCsignalTunedOnData(track);
1835
1836   nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1837
1838   return nSigma;
1839 }
1840
1841 //______________________________________________________________________________
1842 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1843 {
1844   //
1845   // Calculate the number of sigmas in the EMCAL
1846   //
1847
1848   AliVTrack *track=(AliVTrack*)vtrack;
1849
1850   AliVCluster *matchedClus = NULL;
1851
1852   Double_t mom     = -1.;
1853   Double_t pt      = -1.;
1854   Double_t EovP    = -1.;
1855   Double_t fClsE   = -1.;
1856
1857   Int_t nMatchClus = -1;
1858   Int_t charge     = 0;
1859
1860   // Track matching
1861   nMatchClus = track->GetEMCALcluster();
1862   if(nMatchClus > -1){
1863
1864     mom    = track->P();
1865     pt     = track->Pt();
1866     charge = track->Charge();
1867
1868     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1869
1870     if(matchedClus){
1871
1872       // matched cluster is EMCAL
1873       if(matchedClus->IsEMCAL()){
1874
1875         fClsE       = matchedClus->E();
1876         EovP        = fClsE/mom;
1877
1878
1879         // NSigma value really meaningful only for electrons!
1880         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1881       }
1882     }
1883   }
1884
1885   return -999;
1886
1887 }
1888
1889
1890 //______________________________________________________________________________
1891 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1892 {
1893   //
1894   // Compute PID response of 'detCode'
1895   //
1896
1897   switch (detCode){
1898     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1899     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1900     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1901     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1902     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1903     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1904     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1905     default: return kDetNoSignal;
1906   }
1907 }
1908
1909 //______________________________________________________________________________
1910 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1911 {
1912   //
1913   // Compute PID response for the ITS
1914   //
1915
1916   // set flat distribution (no decision)
1917   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1918
1919   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
1920     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
1921
1922   //check for ITS standalone tracks
1923   Bool_t isSA=kTRUE;
1924   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1925
1926   Double_t mom=track->P();
1927   Double_t dedx=track->GetITSsignal();
1928   Double_t momITS=mom;
1929   UChar_t clumap=track->GetITSClusterMap();
1930   Int_t nPointsForPid=0;
1931   for(Int_t i=2; i<6; i++){
1932     if(clumap&(1<<i)) ++nPointsForPid;
1933   }
1934
1935   if(nPointsForPid<3) { // track not to be used for combined PID purposes
1936     //       track->ResetStatus(AliVTrack::kITSpid);
1937     return kDetNoSignal;
1938   }
1939
1940   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1941   for (Int_t j=0; j<nSpecies; j++) {
1942     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1943     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1944     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1945     //TODO: in case of the electron, use the SA parametrisation,
1946     //      this needs to be changed if ITS provides a parametrisation
1947     //      for electrons also for ITS+TPC tracks
1948     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1949     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1950       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1951     } else {
1952       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1953       mismatch=kFALSE;
1954     }
1955
1956     // Check for particles heavier than (AliPID::kSPECIES - 1)
1957     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
1958
1959   }
1960
1961   if (mismatch){
1962     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1963     return kDetNoSignal;
1964   }
1965
1966
1967   return kDetPidOk;
1968 }
1969 //______________________________________________________________________________
1970 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1971 {
1972   //
1973   // Compute PID response for the TPC
1974   //
1975
1976   // set flat distribution (no decision)
1977   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1978
1979   // check quality of the track
1980   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
1981
1982   Double_t dedx=track->GetTPCsignal();
1983   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1984
1985   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1986
1987   Double_t bethe = 0.;
1988   Double_t sigma = 0.;
1989
1990   for (Int_t j=0; j<nSpecies; j++) {
1991     AliPID::EParticleType type=AliPID::EParticleType(j);
1992
1993     bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1994     sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1995
1996     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1997       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1998     } else {
1999       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
2000       mismatch=kFALSE;
2001     }
2002   }
2003
2004   if (mismatch){
2005     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2006     return kDetNoSignal;
2007   }
2008
2009   return kDetPidOk;
2010 }
2011 //______________________________________________________________________________
2012 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2013 {
2014   //
2015   // Compute PID probabilities for TOF
2016   //
2017
2018   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2019
2020   // set flat distribution (no decision)
2021   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2022
2023   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2024   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2025
2026   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2027   for (Int_t j=0; j<nSpecies; j++) {
2028     AliPID::EParticleType type=AliPID::EParticleType(j);
2029     Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor;
2030
2031     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
2032     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
2033     if (TMath::Abs(nsigmas) > (fRange+2)) {
2034       if(nsigmas < fTOFtail)
2035         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
2036       else
2037         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
2038     } else{
2039       if(nsigmas < fTOFtail)
2040         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
2041       else
2042         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
2043     }
2044
2045     if (TMath::Abs(nsigmas)<5.){
2046       Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2047       if (TMath::Abs(nsigmasTPC)>998) mismatch=kFALSE; // if TPC not available we can't check mismatch
2048       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2049     }
2050   }
2051
2052   if (mismatch){
2053     return kDetMismatch;
2054   }
2055
2056   // TODO: Light nuclei
2057
2058   return kDetPidOk;
2059 }
2060 //______________________________________________________________________________
2061 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
2062 {
2063   //
2064   // Compute PID response for the
2065   //
2066
2067   UInt_t TRDslicesForPID[2];
2068   SetTRDSlices(TRDslicesForPID,PIDmethod);
2069   // set flat distribution (no decision)
2070   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2071   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2072
2073   Float_t mom[6]={0.};
2074   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
2075   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2076   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
2077   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2078     mom[ilayer] = track->GetTRDmomentum(ilayer);
2079     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2080       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2081     }
2082   }
2083   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2084   return kDetPidOk;
2085
2086 }
2087 //______________________________________________________________________________
2088 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2089 {
2090   //
2091   // Compute PID response for the EMCAL
2092   //
2093
2094   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2095
2096   AliVCluster *matchedClus = NULL;
2097
2098   Double_t mom     = -1.;
2099   Double_t pt      = -1.;
2100   Double_t EovP    = -1.;
2101   Double_t fClsE   = -1.;
2102
2103   Int_t nMatchClus = -1;
2104   Int_t charge     = 0;
2105
2106   // Track matching
2107   nMatchClus = track->GetEMCALcluster();
2108
2109   if(nMatchClus > -1){
2110
2111     mom    = track->P();
2112     pt     = track->Pt();
2113     charge = track->Charge();
2114
2115     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2116
2117     if(matchedClus){
2118
2119       // matched cluster is EMCAL
2120       if(matchedClus->IsEMCAL()){
2121
2122         fClsE       = matchedClus->E();
2123         EovP        = fClsE/mom;
2124
2125
2126         // compute the probabilities
2127         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
2128
2129           // in case everything is OK
2130           return kDetPidOk;
2131         }
2132       }
2133     }
2134   }
2135
2136   // in all other cases set flat distribution (no decision)
2137   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
2138   return kDetNoSignal;
2139
2140 }
2141 //______________________________________________________________________________
2142 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2143 {
2144   //
2145   // Compute PID response for the PHOS
2146   //
2147
2148   // set flat distribution (no decision)
2149   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2150   return kDetNoSignal;
2151 }
2152 //______________________________________________________________________________
2153 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2154 {
2155   //
2156   // Compute PID response for the HMPID
2157   //
2158   // set flat distribution (no decision)
2159   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2160   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
2161
2162   track->GetHMPIDpid(p);
2163
2164   return kDetPidOk;
2165 }