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(".*(LHC12a17[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="";*/ }
 719
 720   //exception new pp MC productions from 2011
 721   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 722   // exception for 11f1
 723   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 724 }
 725
 726 //______________________________________________________________________________
 727 void AliPIDResponse::SetITSParametrisation()
 728 {
 729   //
 730   // Set the ITS parametrisation
 731   //
 732 }
 733
 734
 735 //______________________________________________________________________________
 736 void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
 737 {
 738   if (h->GetBinContent(binX, binY) <= 1e-4)
 739     return; // Reject bins without content (within some numerical precision) or with strange content
 740
 741   Double_t coord[2] = {0, 0};
 742   coord[0] = h->GetXaxis()->GetBinCenter(binX);
 743   coord[1] = h->GetYaxis()->GetBinCenter(binY);
 744   Double_t binError = h->GetBinError(binX, binY);
 745   if (binError <= 0) {
 746     binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
 747     printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
 748   }
 749   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
 750 }
 751
 752
 753 //______________________________________________________________________________
 754 TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
 755 {
 756   if (!h)
 757     return 0x0;
 758
 759   // Interpolate to finer map
 760   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
 761
 762   Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
 763   Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
 764   Int_t nBinsX = 20;
 765   // 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,
 766   // scale the number of bins correspondingly
 767   Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - lowerMapBoundY) * 40);
 768   Int_t nBinsXrefined = nBinsX * refineFactorX;
 769   Int_t nBinsYrefined = nBinsY * refineFactorY;
 770
 771   TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()),  Form("%s (refined)", h->GetTitle()),
 772                             nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
 773                             nBinsYrefined, lowerMapBoundY, upperMapBoundY);
 774
 775   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
 776     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
 777
 778       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
 779
 780       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
 781       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
 782
 783       /*TODO NOW NOW
 784       linExtrapolation->ClearPoints();
 785
 786       // For interpolation: Just take the corresponding bin from the old histo.
 787       // For extrapolation: take the last available bin from the old histo.
 788       // If the boundaries are to be skipped, also skip the corresponding bins
 789       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
 790       if (oldBinX < 1)
 791         oldBinX = 1;
 792       if (oldBinX > nBinsX)
 793         oldBinX = nBinsX;
 794
 795       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
 796       if (oldBinY < 1)
 797         oldBinY = 1;
 798       if (oldBinY > nBinsY)
 799         oldBinY = nBinsY;
 800
 801       // Neighbours left column
 802       if (oldBinX >= 2) {
 803         if (oldBinY >= 2) {
 804           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
 805         }
 806
 807         AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
 808
 809         if (oldBinY < nBinsY) {
 810           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
 811         }
 812       }
 813
 814       // Neighbours (and point itself) same column
 815       if (oldBinY >= 2) {
 816         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
 817       }
 818
 819       AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
 820
 821       if (oldBinY < nBinsY) {
 822         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
 823       }
 824
 825       // Neighbours right column
 826       if (oldBinX < nBinsX) {
 827         if (oldBinY >= 2) {
 828           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
 829         }
 830
 831         AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
 832
 833         if (oldBinY < nBinsY) {
 834           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
 835         }
 836       }
 837
 838
 839       // Fit 2D-hyperplane
 840       if (linExtrapolation->GetNpoints() <= 0)
 841         continue;
 842
 843       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
 844         continue;
 845
 846       // Fill the bin of the refined histogram with the extrapolated value
 847       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
 848                                  + linExtrapolation->GetParameter(2) * centerY;
 849       */
 850       Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;//TODO NOW NOW NOW
 851       hRefined->SetBinContent(binX, binY, interpolatedValue);
 852     }
 853   }
 854
 855   delete linExtrapolation;
 856
 857   return hRefined;
 858 }
 859
 860 //______________________________________________________________________________
 861 void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
 862                                    Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
 863 {
 864   //
 865   // Load the TPC eta correction maps from the OADB
 866   //
 867
 868   TString dataType = "DATA";
 869   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
 870
 871   if (fIsMC)  {
 872     dataType = "MC";
 873     fRecoPass = 1;
 874
 875     if (fMCperiodTPC.IsNull()) {
 876       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
 877       return;
 878     }
 879
 880     period = fMCperiodTPC;
 881   }
 882
 883   TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), fRecoPass);
 884
 885   AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), fRecoPass));
 886   if (fUseTPCEtaCorrection == kFALSE) {
 887     // Disable eta correction via setting no maps
 888     if (!fTPCResponse.SetEtaCorrMap(0x0))
 889       AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
 890     else
 891       AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
 892
 893     if (!fTPCResponse.SetSigmaParams(0x0, 0))
 894       AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
 895     else
 896       AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
 897
 898     return;
 899   }
 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(), fRecoPass));
 907
 908   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
 909                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
 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) {
 919       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
 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(), fRecoPass));
 957
 958   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
 959                                              Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
 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) {
 969       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
 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   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1100   Bool_t found=kFALSE;
1101   //
1102   //set the new PID splines
1103   //
1104   if (fArrPidResponseMaster){
1105     Int_t recopass = fRecoPass;
1106     if(fTuneMConData) recopass = fRecoPassUser;
1107     //for MC don't use period information
1108     //if (fIsMC) period="[A-Z0-9]*";
1109     //for MC use MC period information
1110     //pattern for the default entry (valid for all particles)
1111     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1112
1113     //find particle id ang gain scenario
1114     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1115     {
1116       TObject *grAll=NULL;
1117       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1118       gainScenario.ToUpper();
1119       //loop over entries and filter them
1120       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1121       {
1122         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1123         if (responseFunction==NULL) continue;
1124         TString responseName=responseFunction->GetName();
1125
1126         if (!reg.MatchB(responseName)) continue;
1127
1128         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1129         TObject* tmp=NULL;
1130         tmp=arr->At(1); if (!tmp) continue;
1131         TString particleName=tmp->GetName();
1132         tmp=arr->At(3); if (!tmp) continue;
1133         TString gainScenarioName=tmp->GetName();
1134         delete arr;
1135         if (particleName.IsNull()) continue;
1136         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1137         else
1138         {
1139           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1140           {
1141             TString particle=AliPID::ParticleName(ispec);
1142             particle.ToUpper();
1143             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1144             if ( particle == particleName && gainScenario == gainScenarioName )
1145             {
1146               fTPCResponse.SetResponseFunction( responseFunction,
1147                                                 (AliPID::EParticleType)ispec,
1148                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1149               fTPCResponse.SetUseDatabase(kTRUE);
1150               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
1151               found=kTRUE;
1152               // overwrite default with proton spline (for light nuclei)
1153               if (ispec==AliPID::kProton) grAll=responseFunction;
1154               break;
1155             }
1156           }
1157         }
1158       }
1159       if (grAll)
1160       {
1161         for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1162         {
1163           if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1164                                                  (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1165           {
1166               fTPCResponse.SetResponseFunction( grAll,
1167                                                 (AliPID::EParticleType)ispec,
1168                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1169               fTPCResponse.SetUseDatabase(kTRUE);
1170               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1171               found=kTRUE;
1172           }
1173         }
1174       }
1175     }
1176   }
1177   else AliInfo("no fArrPidResponseMaster");
1178
1179   if (!found){
1180     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1181   }
1182
1183   //
1184   // Setup resolution parametrisation
1185   //
1186
1187   //default
1188   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1189
1190   if (fRun>=122195){
1191     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1192   }
1193
1194   if (fRun>=186636){
1195 //   if (fRun>=188356){
1196     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1197   }
1198
1199   if (fArrPidResponseMaster)
1200   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1201
1202   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1203
1204   //read in the voltage map
1205   TVectorF* gsm = 0x0;
1206   if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1207   if (gsm)
1208   {
1209     fTPCResponse.SetVoltageMap(*gsm);
1210     TString vals;
1211     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1212     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1213     AliInfo(vals.Data());
1214     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1215     AliInfo(vals.Data());
1216     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1217     AliInfo(vals.Data());
1218     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1219     AliInfo(vals.Data());
1220   }
1221   else AliInfo("no voltage map, ideal default assumed");
1222 }
1223
1224 //______________________________________________________________________________
1225 void AliPIDResponse::SetTRDPidResponseMaster()
1226 {
1227   //
1228   // Load the TRD pid params and references from the OADB
1229   //
1230   if(fTRDPIDResponseObject) return;
1231   AliOADBContainer contParams("contParams");
1232
1233   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1234   if(statusResponse){
1235     AliError("Failed initializing PID Response Object from OADB");
1236   } else {
1237     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1238     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1239     if(!fTRDPIDResponseObject){
1240       AliError(Form("TRD Response not found in run %d", fRun));
1241     }
1242   }
1243 }
1244
1245 //______________________________________________________________________________
1246 void AliPIDResponse::InitializeTRDResponse(){
1247   //
1248   // Set PID Params and references to the TRD PID response
1249   //
1250     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1251 }
1252
1253 //______________________________________________________________________________
1254 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1255
1256     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1257         // backward compatibility for setting with 8 slices
1258         TRDslicesForPID[0] = 0;
1259         TRDslicesForPID[1] = 7;
1260     }
1261     else{
1262         if(method==AliTRDPIDResponse::kLQ1D){
1263             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1264             TRDslicesForPID[1] = 0;
1265         }
1266         if(method==AliTRDPIDResponse::kLQ2D){
1267             TRDslicesForPID[0] = 1;
1268             TRDslicesForPID[1] = 7;
1269         }
1270     }
1271     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1272 }
1273
1274 //______________________________________________________________________________
1275 void AliPIDResponse::SetTOFPidResponseMaster()
1276 {
1277   //
1278   // Load the TOF pid params from the OADB
1279   //
1280
1281   if (fTOFPIDParams) delete fTOFPIDParams;
1282   fTOFPIDParams=NULL;
1283
1284   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1285   if (oadbf && oadbf->IsOpen()) {
1286     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1287     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1288     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1289     oadbf->Close();
1290     delete oadbc;
1291   }
1292   delete oadbf;
1293
1294   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1295 }
1296
1297 //______________________________________________________________________________
1298 void AliPIDResponse::InitializeTOFResponse(){
1299   //
1300   // Set PID Params to the TOF PID response
1301   //
1302
1303   AliInfo("TOF PID Params loaded from OADB");
1304   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1305   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1306   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1307                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1308
1309   for (Int_t i=0;i<4;i++) {
1310     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1311   }
1312   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1313
1314   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1315   Float_t t0Spread[4];
1316   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1317   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]));
1318   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1319   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1320   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1321     fResT0AC=t0Spread[3];
1322     fResT0A=TMath::Sqrt(a);
1323     fResT0C=TMath::Sqrt(c);
1324   } else {
1325     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1326     fResT0A=75.;
1327     fResT0C=65.;
1328     fResT0AC=55.;
1329   }
1330   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1331
1332 }
1333
1334
1335 //______________________________________________________________________________
1336 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1337   //
1338   // Check whether track is identified as electron under a given electron efficiency hypothesis
1339     //
1340
1341   Double_t probs[AliPID::kSPECIES];
1342   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1343
1344   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1345   // Take mean of the TRD momenta in the given tracklets
1346   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1347   Int_t nmomenta = 0;
1348   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1349     if(vtrack->GetTRDmomentum(iPl) > 0.){
1350       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1351     }
1352   }
1353   p = TMath::Mean(nmomenta, trdmomenta);
1354
1355   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1356 }
1357
1358 //______________________________________________________________________________
1359 void AliPIDResponse::SetEMCALPidResponseMaster()
1360 {
1361   //
1362   // Load the EMCAL pid response functions from the OADB
1363   //
1364   TObjArray* fEMCALPIDParamsRun      = NULL;
1365   TObjArray* fEMCALPIDParamsPass     = NULL;
1366
1367   if(fEMCALPIDParams) return;
1368   AliOADBContainer contParams("contParams");
1369
1370   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1371   if(statusPars){
1372     AliError("Failed initializing PID Params from OADB");
1373   }
1374   else {
1375     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1376
1377     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1378     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1379     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1380
1381     if(!fEMCALPIDParams){
1382       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1383       AliInfo("Will take the standard LHC11d instead ...");
1384
1385       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1386       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1387       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1388
1389       if(!fEMCALPIDParams){
1390         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1391       }
1392     }
1393   }
1394 }
1395
1396 //______________________________________________________________________________
1397 void AliPIDResponse::InitializeEMCALResponse(){
1398   //
1399   // Set PID Params to the EMCAL PID response
1400   //
1401   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1402
1403 }
1404
1405 //______________________________________________________________________________
1406 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1407 {
1408   //
1409   // create detector PID information and setup the transient pointer in the track
1410   //
1411
1412   // check if detector number is inside accepted range
1413   if (detector == kNdetectors) return;
1414
1415   // get detector pid
1416   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1417   if (!detPID) {
1418     detPID=new AliDetectorPID;
1419     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1420   }
1421
1422   //check if values exist
1423   if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return;
1424
1425   //TODO: which particles to include? See also the loops below...
1426   Double_t values[AliPID::kSPECIESC]={0};
1427
1428   //nsigmas
1429   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1430     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1431   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC);
1432
1433   //probabilities
1434   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1435   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1436 }
1437
1438 //______________________________________________________________________________
1439 void AliPIDResponse::FillTrackDetectorPID()
1440 {
1441   //
1442   // create detector PID information and setup the transient pointer in the track
1443   //
1444
1445   if (!fCurrentEvent) return;
1446
1447   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1448     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1449     if (!track) continue;
1450
1451     for (Int_t idet=0; idet<kNdetectors; ++idet){
1452       FillTrackDetectorPID(track, (EDetector)idet);
1453     }
1454   }
1455 }
1456
1457 //______________________________________________________________________________
1458 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1459   //
1460   // Set TOF response function
1461   // Input option for event_time used
1462   //
1463
1464     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1465     if(t0spread < 10) t0spread = 80;
1466
1467     // T0 from TOF algorithm
1468
1469     Bool_t flagT0TOF=kFALSE;
1470     Bool_t flagT0T0=kFALSE;
1471     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1472     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1473     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1474
1475     // T0-TOF arrays
1476     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1477     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1478     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1479       estimatedT0event[i]=0.0;
1480       estimatedT0resolution[i]=0.0;
1481       startTimeMask[i] = 0;
1482     }
1483
1484     Float_t resT0A=fResT0A;
1485     Float_t resT0C=fResT0C;
1486     Float_t resT0AC=fResT0AC;
1487     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1488         flagT0T0=kTRUE;
1489     }
1490
1491
1492     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1493
1494     if (tofHeader) { // read global info and T0-TOF
1495       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1496       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1497       if(t0spread < 10) t0spread = 80;
1498
1499       flagT0TOF=kTRUE;
1500       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1501         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1502         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1503         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1504       }
1505
1506       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1507       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1508       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1509       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1510         Int_t icurrent = (Int_t)ibin->GetAt(j);
1511         startTime[icurrent]=t0Bin->GetAt(j);
1512         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1513         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1514       }
1515     }
1516
1517     // for cut of 3 sigma on t0 spread
1518     Float_t t0cut = 3 * t0spread;
1519     if(t0cut < 500) t0cut = 500;
1520
1521     if(option == kFILL_T0){ // T0-FILL is used
1522         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1523           estimatedT0event[i]=0.0;
1524           estimatedT0resolution[i]=t0spread;
1525         }
1526         fTOFResponse.SetT0event(estimatedT0event);
1527         fTOFResponse.SetT0resolution(estimatedT0resolution);
1528     }
1529
1530     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1531         if(flagT0TOF){
1532             fTOFResponse.SetT0event(startTime);
1533             fTOFResponse.SetT0resolution(startTimeRes);
1534             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1535               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1536               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1537             }
1538         }
1539         else{
1540             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1541               estimatedT0event[i]=0.0;
1542               estimatedT0resolution[i]=t0spread;
1543               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1544             }
1545             fTOFResponse.SetT0event(estimatedT0event);
1546             fTOFResponse.SetT0resolution(estimatedT0resolution);
1547         }
1548     }
1549     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1550         Float_t t0AC=-10000;
1551         Float_t t0A=-10000;
1552         Float_t t0C=-10000;
1553         if(flagT0T0){
1554             t0A= vevent->GetT0TOF()[1];
1555             t0C= vevent->GetT0TOF()[2];
1556         //      t0AC= vevent->GetT0TOF()[0];
1557         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1558         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1559         t0AC /= resT0AC*resT0AC;
1560         }
1561
1562         Float_t t0t0Best = 0;
1563         Float_t t0t0BestRes = 9999;
1564         Int_t t0used=0;
1565         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1566             t0t0Best = t0AC;
1567             t0t0BestRes = resT0AC;
1568             t0used=6;
1569         }
1570         else if(TMath::Abs(t0C) < t0cut){
1571             t0t0Best = t0C;
1572             t0t0BestRes = resT0C;
1573             t0used=4;
1574         }
1575         else if(TMath::Abs(t0A) < t0cut){
1576             t0t0Best = t0A;
1577             t0t0BestRes = resT0A;
1578             t0used=2;
1579         }
1580
1581         if(flagT0TOF){ // if T0-TOF info is available
1582             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1583                 if(t0t0BestRes < 999){
1584                   if(startTimeRes[i] < t0spread){
1585                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1586                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1587                     estimatedT0event[i]=t0best / wtot;
1588                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1589                     startTimeMask[i] = t0used+1;
1590                   }
1591                   else {
1592                     estimatedT0event[i]=t0t0Best;
1593                     estimatedT0resolution[i]=t0t0BestRes;
1594                     startTimeMask[i] = t0used;
1595                   }
1596                 }
1597                 else{
1598                   estimatedT0event[i]=startTime[i];
1599                   estimatedT0resolution[i]=startTimeRes[i];
1600                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1601                 }
1602                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1603             }
1604             fTOFResponse.SetT0event(estimatedT0event);
1605             fTOFResponse.SetT0resolution(estimatedT0resolution);
1606         }
1607         else{ // if no T0-TOF info is available
1608             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1609               fTOFResponse.SetT0binMask(i,t0used);
1610               if(t0t0BestRes < 999){
1611                 estimatedT0event[i]=t0t0Best;
1612                 estimatedT0resolution[i]=t0t0BestRes;
1613               }
1614               else{
1615                 estimatedT0event[i]=0.0;
1616                 estimatedT0resolution[i]=t0spread;
1617               }
1618             }
1619             fTOFResponse.SetT0event(estimatedT0event);
1620             fTOFResponse.SetT0resolution(estimatedT0resolution);
1621         }
1622     }
1623
1624     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1625         Float_t t0AC=-10000;
1626         Float_t t0A=-10000;
1627         Float_t t0C=-10000;
1628         if(flagT0T0){
1629             t0A= vevent->GetT0TOF()[1];
1630             t0C= vevent->GetT0TOF()[2];
1631         //      t0AC= vevent->GetT0TOF()[0];
1632         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1633         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1634         t0AC /= resT0AC*resT0AC;
1635         }
1636
1637         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1638             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1639               estimatedT0event[i]=t0AC;
1640               estimatedT0resolution[i]=resT0AC;
1641               fTOFResponse.SetT0binMask(i,6);
1642             }
1643         }
1644         else if(TMath::Abs(t0C) < t0cut){
1645             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1646               estimatedT0event[i]=t0C;
1647               estimatedT0resolution[i]=resT0C;
1648               fTOFResponse.SetT0binMask(i,4);
1649             }
1650         }
1651         else if(TMath::Abs(t0A) < t0cut){
1652             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1653               estimatedT0event[i]=t0A;
1654               estimatedT0resolution[i]=resT0A;
1655               fTOFResponse.SetT0binMask(i,2);
1656             }
1657         }
1658         else{
1659             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1660               estimatedT0event[i]=0.0;
1661               estimatedT0resolution[i]=t0spread;
1662               fTOFResponse.SetT0binMask(i,0);
1663             }
1664         }
1665         fTOFResponse.SetT0event(estimatedT0event);
1666         fTOFResponse.SetT0resolution(estimatedT0resolution);
1667     }
1668     delete [] startTime;
1669     delete [] startTimeRes;
1670     delete [] startTimeMask;
1671     delete [] estimatedT0event;
1672     delete [] estimatedT0resolution;
1673 }
1674
1675 //______________________________________________________________________________
1676 // private non cached versions of the PID calculation
1677 //
1678
1679
1680 //______________________________________________________________________________
1681 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
1682 {
1683   //
1684   // NumberOfSigmas for 'detCode'
1685   //
1686
1687   switch (detCode){
1688     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1689     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1690     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1691     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1692     default: return -999.;
1693   }
1694
1695 }
1696
1697
1698
1699 //______________________________________________________________________________
1700 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1701 {
1702   //
1703   // Calculate the number of sigmas in the ITS
1704   //
1705
1706   AliVTrack *track=(AliVTrack*)vtrack;
1707
1708   Float_t dEdx=track->GetITSsignal();
1709   if (dEdx<=0) return -999.;
1710
1711   UChar_t clumap=track->GetITSClusterMap();
1712   Int_t nPointsForPid=0;
1713   for(Int_t i=2; i<6; i++){
1714     if(clumap&(1<<i)) ++nPointsForPid;
1715   }
1716   Float_t mom=track->P();
1717
1718   //check for ITS standalone tracks
1719   Bool_t isSA=kTRUE;
1720   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1721
1722   //TODO: in case of the electron, use the SA parametrisation,
1723   //      this needs to be changed if ITS provides a parametrisation
1724   //      for electrons also for ITS+TPC tracks
1725   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1726 }
1727
1728 //______________________________________________________________________________
1729 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1730 {
1731   //
1732   // Calculate the number of sigmas in the TPC
1733   //
1734
1735   AliVTrack *track=(AliVTrack*)vtrack;
1736
1737   Double_t nSigma = -999.;
1738
1739   //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
1740   if (fTuneMConData) {
1741     Double_t mom  = track->GetTPCmomentum();
1742     Double_t sig  = this->GetTPCsignalTunedOnData(track);
1743     UInt_t   sigN = track->GetTPCsignalN();
1744
1745     if (sigN>0)
1746       nSigma = fTPCResponse.GetNumberOfSigmas(mom, sig, sigN, type);
1747   }
1748   else {
1749     nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1750   }
1751
1752   return nSigma;
1753 }
1754
1755 //______________________________________________________________________________
1756 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1757 {
1758   //
1759   // Calculate the number of sigmas in the EMCAL
1760   //
1761
1762   AliVTrack *track=(AliVTrack*)vtrack;
1763
1764   AliVCluster *matchedClus = NULL;
1765
1766   Double_t mom     = -1.;
1767   Double_t pt      = -1.;
1768   Double_t EovP    = -1.;
1769   Double_t fClsE   = -1.;
1770
1771   Int_t nMatchClus = -1;
1772   Int_t charge     = 0;
1773
1774   // Track matching
1775   nMatchClus = track->GetEMCALcluster();
1776   if(nMatchClus > -1){
1777
1778     mom    = track->P();
1779     pt     = track->Pt();
1780     charge = track->Charge();
1781
1782     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1783
1784     if(matchedClus){
1785
1786       // matched cluster is EMCAL
1787       if(matchedClus->IsEMCAL()){
1788
1789         fClsE       = matchedClus->E();
1790         EovP        = fClsE/mom;
1791
1792
1793         // NSigma value really meaningful only for electrons!
1794         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1795       }
1796     }
1797   }
1798
1799   return -999;
1800
1801 }
1802
1803
1804 //______________________________________________________________________________
1805 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1806 {
1807   //
1808   // Compute PID response of 'detCode'
1809   //
1810
1811   switch (detCode){
1812     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1813     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1814     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1815     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1816     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1817     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1818     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1819     default: return kDetNoSignal;
1820   }
1821 }
1822
1823 //______________________________________________________________________________
1824 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1825 {
1826   //
1827   // Compute PID response for the ITS
1828   //
1829
1830   // look for cached value first
1831   // only the non SA tracks are cached
1832   if (track->GetDetectorPID()){
1833     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1834   }
1835
1836   // set flat distribution (no decision)
1837   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1838
1839   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
1840     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
1841
1842   //check for ITS standalone tracks
1843   Bool_t isSA=kTRUE;
1844   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1845
1846   Double_t mom=track->P();
1847   Double_t dedx=track->GetITSsignal();
1848   Double_t momITS=mom;
1849   UChar_t clumap=track->GetITSClusterMap();
1850   Int_t nPointsForPid=0;
1851   for(Int_t i=2; i<6; i++){
1852     if(clumap&(1<<i)) ++nPointsForPid;
1853   }
1854
1855   if(nPointsForPid<3) { // track not to be used for combined PID purposes
1856     //       track->ResetStatus(AliVTrack::kITSpid);
1857     return kDetNoSignal;
1858   }
1859
1860   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1861   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
1862     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1863     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1864     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1865     //TODO: in case of the electron, use the SA parametrisation,
1866     //      this needs to be changed if ITS provides a parametrisation
1867     //      for electrons also for ITS+TPC tracks
1868     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1869     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1870       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1871     } else {
1872       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1873       mismatch=kFALSE;
1874     }
1875
1876     // Check for particles heavier than (AliPID::kSPECIES - 1)
1877     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
1878
1879   }
1880
1881   if (mismatch){
1882     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1883     return kDetNoSignal;
1884   }
1885
1886
1887   return kDetPidOk;
1888 }
1889 //______________________________________________________________________________
1890 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1891 {
1892   //
1893   // Compute PID response for the TPC
1894   //
1895
1896   // set flat distribution (no decision)
1897   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1898
1899   // check quality of the track
1900   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
1901
1902   Double_t dedx=track->GetTPCsignal();
1903   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1904
1905   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1906
1907   Double_t bethe = 0.;
1908   Double_t sigma = 0.;
1909
1910   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1911     AliPID::EParticleType type=AliPID::EParticleType(j);
1912
1913     //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
1914     if (fTuneMConData) {
1915       Double_t mom = track->GetTPCmomentum();
1916       bethe=fTPCResponse.GetExpectedSignal(mom,type);
1917       sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
1918     }
1919     else {
1920       bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1921       sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1922     }
1923     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1924       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1925     } else {
1926       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1927       mismatch=kFALSE;
1928     }
1929   }
1930
1931   if (mismatch){
1932     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1933     return kDetNoSignal;
1934   }
1935
1936   return kDetPidOk;
1937 }
1938 //______________________________________________________________________________
1939 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1940 {
1941   //
1942   // Compute PID response for the
1943   //
1944
1945   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1946
1947   // set flat distribution (no decision)
1948   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1949
1950   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
1951   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
1952
1953   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
1954   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1955     AliPID::EParticleType type=AliPID::EParticleType(j);
1956     Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor;
1957
1958     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1959     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1960     if (TMath::Abs(nsigmas) > (fRange+2)) {
1961       if(nsigmas < fTOFtail)
1962         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1963       else
1964         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1965     } else{
1966       if(nsigmas < fTOFtail)
1967         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1968       else
1969         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1970     }
1971
1972     if (TMath::Abs(nsigmas)<5.){
1973       Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
1974       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
1975     }
1976   }
1977
1978   if (mismatch){
1979     return kDetMismatch;
1980   }
1981
1982   // TODO: Light nuclei
1983
1984   return kDetPidOk;
1985 }
1986 //______________________________________________________________________________
1987 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1988 {
1989   //
1990   // Compute PID response for the
1991   //
1992
1993   UInt_t TRDslicesForPID[2];
1994   SetTRDSlices(TRDslicesForPID,PIDmethod);
1995   // set flat distribution (no decision)
1996   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1997   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
1998
1999   Float_t mom[6]={0.};
2000   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
2001   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2002   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
2003   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2004     mom[ilayer] = track->GetTRDmomentum(ilayer);
2005     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2006       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2007     }
2008   }
2009   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2010   return kDetPidOk;
2011
2012 }
2013 //______________________________________________________________________________
2014 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2015 {
2016   //
2017   // Compute PID response for the EMCAL
2018   //
2019
2020   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2021
2022   AliVCluster *matchedClus = NULL;
2023
2024   Double_t mom     = -1.;
2025   Double_t pt      = -1.;
2026   Double_t EovP    = -1.;
2027   Double_t fClsE   = -1.;
2028
2029   Int_t nMatchClus = -1;
2030   Int_t charge     = 0;
2031
2032   // Track matching
2033   nMatchClus = track->GetEMCALcluster();
2034
2035   if(nMatchClus > -1){
2036
2037     mom    = track->P();
2038     pt     = track->Pt();
2039     charge = track->Charge();
2040
2041     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2042
2043     if(matchedClus){
2044
2045       // matched cluster is EMCAL
2046       if(matchedClus->IsEMCAL()){
2047
2048         fClsE       = matchedClus->E();
2049         EovP        = fClsE/mom;
2050
2051
2052         // compute the probabilities
2053         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
2054
2055           // in case everything is OK
2056           return kDetPidOk;
2057         }
2058       }
2059     }
2060   }
2061
2062   // in all other cases set flat distribution (no decision)
2063   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
2064   return kDetNoSignal;
2065
2066 }
2067 //______________________________________________________________________________
2068 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2069 {
2070   //
2071   // Compute PID response for the PHOS
2072   //
2073
2074   // set flat distribution (no decision)
2075   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2076   return kDetNoSignal;
2077 }
2078 //______________________________________________________________________________
2079 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2080 {
2081   //
2082   // Compute PID response for the HMPID
2083   //
2084   // set flat distribution (no decision)
2085   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2086   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
2087
2088   track->GetHMPIDpid(p);
2089
2090   return kDetPidOk;
2091 }