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