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