STEER/STEERBase/AliPIDResponse.cxx

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