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