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 <TSpline.h>
  31 #include <TFile.h>
  32
  33 #include <AliVEvent.h>
  34 #include <AliVTrack.h>
  35 #include <AliLog.h>
  36 #include <AliPID.h>
  37 #include <AliOADBContainer.h>
  38 #include <AliTRDPIDParams.h>
  39 #include <AliTRDPIDReference.h>
  40
  41 #include "AliPIDResponse.h"
  42
  43 #include "AliCentrality.h"
  44
  45 ClassImp(AliPIDResponse);
  46
  47 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
  48 TNamed("PIDResponse","PIDResponse"),
  49 fITSResponse(isMC),
  50 fTPCResponse(),
  51 fTRDResponse(),
  52 fTOFResponse(),
  53 fEMCALResponse(),
  54 fRange(5.),
  55 fITSPIDmethod(kITSTruncMean),
  56 fIsMC(isMC),
  57 fOADBPath(),
  58 fBeamType("PP"),
  59 fLHCperiod(),
  60 fMCperiodTPC(),
  61 fMCperiodUser(),
  62 fCurrentFile(),
  63 fRecoPass(0),
  64 fRecoPassUser(-1),
  65 fRun(0),
  66 fOldRun(0),
  67 fArrPidResponseMaster(0x0),
  68 fResolutionCorrection(0x0),
  69 fTRDPIDParams(0x0),
  70 fTRDPIDReference(0x0),
  71 fTOFTimeZeroType(kBest_T0),
  72 fTOFres(100.),
  73 fEMCALPIDParams(0x0),
  74 fCurrentEvent(0x0),
  75 fCurrCentrality(0.0)
  76 {
  77   //
  78   // default ctor
  79   //
  80   AliLog::SetClassDebugLevel("AliPIDResponse",10);
  81   AliLog::SetClassDebugLevel("AliESDpid",10);
  82   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
  83
  84   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
  85 }
  86
  87 //______________________________________________________________________________
  88 AliPIDResponse::~AliPIDResponse()
  89 {
  90   //
  91   // dtor
  92   //
  93   delete fArrPidResponseMaster;
  94   delete fTRDPIDParams;
  95   delete fTRDPIDReference;
  96 }
  97
  98 //______________________________________________________________________________
  99 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
 100 TNamed(other),
 101 fITSResponse(other.fITSResponse),
 102 fTPCResponse(other.fTPCResponse),
 103 fTRDResponse(other.fTRDResponse),
 104 fTOFResponse(other.fTOFResponse),
 105 fEMCALResponse(other.fEMCALResponse),
 106 fRange(other.fRange),
 107 fITSPIDmethod(other.fITSPIDmethod),
 108 fIsMC(other.fIsMC),
 109 fOADBPath(other.fOADBPath),
 110 fBeamType("PP"),
 111 fLHCperiod(),
 112 fMCperiodTPC(),
 113 fMCperiodUser(other.fMCperiodUser),
 114 fCurrentFile(),
 115 fRecoPass(0),
 116 fRecoPassUser(other.fRecoPassUser),
 117 fRun(0),
 118 fOldRun(0),
 119 fArrPidResponseMaster(0x0),
 120 fResolutionCorrection(0x0),
 121 fTRDPIDParams(0x0),
 122 fTRDPIDReference(0x0),
 123 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
 124 fTOFres(100.),
 125 fEMCALPIDParams(0x0),
 126 fCurrentEvent(0x0),
 127 fCurrCentrality(0.0)
 128 {
 129   //
 130   // copy ctor
 131   //
 132   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 133 }
 134
 135 //______________________________________________________________________________
 136 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
 137 {
 138   //
 139   // copy ctor
 140   //
 141   if(this!=&other) {
 142     delete fArrPidResponseMaster;
 143     TNamed::operator=(other);
 144     fITSResponse=other.fITSResponse;
 145     fTPCResponse=other.fTPCResponse;
 146     fTRDResponse=other.fTRDResponse;
 147     fTOFResponse=other.fTOFResponse;
 148     fEMCALResponse=other.fEMCALResponse;
 149     fRange=other.fRange;
 150     fITSPIDmethod=other.fITSPIDmethod;
 151     fOADBPath=other.fOADBPath;
 152     fIsMC=other.fIsMC;
 153     fBeamType="PP";
 154     fLHCperiod="";
 155     fMCperiodTPC="";
 156     fMCperiodUser=other.fMCperiodUser;
 157     fCurrentFile="";
 158     fRecoPass=0;
 159     fRecoPassUser=other.fRecoPassUser;
 160     fRun=0;
 161     fOldRun=0;
 162     fArrPidResponseMaster=0x0;
 163     fResolutionCorrection=0x0;
 164     fTRDPIDParams=0x0;
 165     fTRDPIDReference=0x0;
 166     fEMCALPIDParams=0x0;
 167     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 168     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
 169     fTOFres=100.;
 170     fCurrentEvent=other.fCurrentEvent;
 171   }
 172   return *this;
 173 }
 174
 175 //______________________________________________________________________________
 176 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
 177 {
 178   //
 179   // NumberOfSigmas for 'detCode'
 180   //
 181
 182   switch (detCode){
 183     case kDetITS: return NumberOfSigmasITS(track, type); break;
 184     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
 185     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
 186 //     case kDetTRD: return ComputeTRDProbability(track, type); break;
 187 //     case kDetPHOS: return ComputePHOSProbability(track, type); break;
 188 //     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
 189 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
 190     default: return -999.;
 191   }
 192
 193 }
 194
 195 //______________________________________________________________________________
 196 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
 197
 198   AliVCluster *matchedClus = NULL;
 199
 200   Double_t mom     = -1.;
 201   Double_t pt      = -1.;
 202   Double_t EovP    = -1.;
 203   Double_t fClsE   = -1.;
 204
 205   Int_t nMatchClus = -1;
 206   Int_t charge     = 0;
 207
 208   // Track matching
 209   nMatchClus = track->GetEMCALcluster();
 210   if(nMatchClus > -1){
 211
 212     mom    = track->P();
 213     pt     = track->Pt();
 214     charge = track->Charge();
 215
 216     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 217
 218     if(matchedClus){
 219
 220     // matched cluster is EMCAL
 221     if(matchedClus->IsEMCAL()){
 222
 223       fClsE       = matchedClus->E();
 224       EovP        = fClsE/mom;
 225
 226
 227       // NSigma value really meaningful only for electrons!
 228       return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
 229     }
 230   }
 231   }
 232
 233   return -999;
 234
 235 }
 236
 237 //______________________________________________________________________________
 238 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 239 {
 240   //
 241   // Compute PID response of 'detCode'
 242   //
 243
 244   switch (detCode){
 245     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
 246     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
 247     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
 248     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
 249     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
 250     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
 251     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
 252     default: return kDetNoSignal;
 253   }
 254 }
 255
 256 //______________________________________________________________________________
 257 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 258 {
 259   //
 260   // Compute PID response for the ITS
 261   //
 262
 263   // set flat distribution (no decision)
 264   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 265
 266   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
 267     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
 268
 269   Double_t mom=track->P();
 270   Double_t dedx=track->GetITSsignal();
 271   Bool_t isSA=kTRUE;
 272   Double_t momITS=mom;
 273   ULong_t trStatus=track->GetStatus();
 274   if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
 275   UChar_t clumap=track->GetITSClusterMap();
 276   Int_t nPointsForPid=0;
 277   for(Int_t i=2; i<6; i++){
 278     if(clumap&(1<<i)) ++nPointsForPid;
 279   }
 280
 281   if(nPointsForPid<3) { // track not to be used for combined PID purposes
 282     //       track->ResetStatus(AliVTrack::kITSpid);
 283     return kDetNoSignal;
 284   }
 285
 286   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
 287   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 288     Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
 289     Double_t bethe=fITSResponse.Bethe(momITS,mass);
 290     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
 291     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 292       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 293     } else {
 294       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 295       mismatch=kFALSE;
 296     }
 297
 298     // Check for particles heavier than (AliPID::kSPECIES - 1)
 299     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
 300
 301   }
 302
 303   if (mismatch){
 304     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 305     return kDetNoSignal;
 306   }
 307
 308
 309   return kDetPidOk;
 310 }
 311 //______________________________________________________________________________
 312 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 313 {
 314   //
 315   // Compute PID response for the TPC
 316   //
 317
 318   // set flat distribution (no decision)
 319   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 320
 321   // check quality of the track
 322   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
 323
 324   Double_t mom = track->GetTPCmomentum();
 325
 326   Double_t dedx=track->GetTPCsignal();
 327   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
 328
 329   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 330     AliPID::EParticleType type=AliPID::EParticleType(j);
 331     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 332     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
 333     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 334       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 335     } else {
 336       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 337       mismatch=kFALSE;
 338     }
 339
 340     // TODO: Light nuclei, also in TPC pid response
 341
 342     // Check for particles heavier than (AliPID::kSPECIES - 1)
 343 //     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
 344
 345   }
 346
 347   if (mismatch){
 348     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 349     return kDetNoSignal;
 350   }
 351
 352   return kDetPidOk;
 353 }
 354 //______________________________________________________________________________
 355 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 356 {
 357   //
 358   // Compute PID response for the
 359   //
 360
 361   // set flat distribution (no decision)
 362   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 363
 364   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
 365   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
 366
 367   Double_t time[AliPID::kSPECIESN];
 368   track->GetIntegratedTimes(time);
 369
 370   Double_t sigma[AliPID::kSPECIES];
 371   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
 372     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
 373   }
 374
 375   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
 376   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 377     AliPID::EParticleType type=AliPID::EParticleType(j);
 378     Double_t nsigmas=NumberOfSigmasTOF(track,type);
 379
 380     Double_t sig = sigma[j];
 381     if (TMath::Abs(nsigmas) > (fRange+2)) {
 382       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
 383     } else
 384       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
 385
 386     if (TMath::Abs(nsigmas)<5.){
 387       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
 388       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
 389     }
 390   }
 391
 392   if (mismatch){
 393     return kDetMismatch;
 394   }
 395
 396     // TODO: Light nuclei
 397
 398   return kDetPidOk;
 399 }
 400 //______________________________________________________________________________
 401 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 402 {
 403   //
 404   // Compute PID response for the
 405   //
 406
 407   // set flat distribution (no decision)
 408   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 409   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
 410
 411   Float_t mom[6];
 412   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
 413   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
 414   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
 415   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
 416     mom[ilayer] = track->GetTRDmomentum(ilayer);
 417     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
 418       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
 419     }
 420   }
 421   fTRDResponse.GetResponse(nslices, dedx, mom, p);
 422   return kDetPidOk;
 423 }
 424 //______________________________________________________________________________
 425 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 426 {
 427   //
 428   // Compute PID response for the EMCAL
 429   //
 430
 431   AliVCluster *matchedClus = NULL;
 432
 433   Double_t mom     = -1.;
 434   Double_t pt      = -1.;
 435   Double_t EovP    = -1.;
 436   Double_t fClsE   = -1.;
 437
 438   Int_t nMatchClus = -1;
 439   Int_t charge     = 0;
 440
 441   // Track matching
 442   nMatchClus = track->GetEMCALcluster();
 443
 444   if(nMatchClus > -1){
 445
 446     mom    = track->P();
 447     pt     = track->Pt();
 448     charge = track->Charge();
 449
 450     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 451
 452     if(matchedClus){
 453
 454     // matched cluster is EMCAL
 455     if(matchedClus->IsEMCAL()){
 456
 457       fClsE       = matchedClus->E();
 458       EovP        = fClsE/mom;
 459
 460
 461       // compute the probabilities
 462       if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
 463
 464         // in case everything is OK
 465         return kDetPidOk;
 466
 467       }
 468     }
 469   }
 470   }
 471
 472   // in all other cases set flat distribution (no decision)
 473   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 474   return kDetNoSignal;
 475
 476 }
 477 //______________________________________________________________________________
 478 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
 479 {
 480   //
 481   // Compute PID response for the PHOS
 482   //
 483
 484   // set flat distribution (no decision)
 485   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 486   return kDetNoSignal;
 487 }
 488 //______________________________________________________________________________
 489 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 490 {
 491   //
 492   // Compute PID response for the HMPID
 493   //
 494
 495   // set flat distribution (no decision)
 496   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 497   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
 498
 499   track->GetHMPIDpid(p);
 500
 501   return kDetPidOk;
 502 }
 503
 504 //______________________________________________________________________________
 505 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
 506 {
 507   //
 508   // Apply settings for the current event
 509   //
 510   fRecoPass=pass;
 511
 512   fCurrentEvent=0x0;
 513   if (!event) return;
 514   fCurrentEvent=event;
 515   fRun=event->GetRunNumber();
 516
 517   if (fRun!=fOldRun){
 518     ExecNewRun();
 519     fOldRun=fRun;
 520   }
 521
 522   //TPC resolution parametrisation PbPb
 523   if ( fResolutionCorrection ){
 524     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
 525     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
 526   }
 527
 528   //TOF resolution
 529   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
 530
 531   // Get and set centrality
 532   AliCentrality *centrality = event->GetCentrality();
 533   if(centrality){
 534     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
 535   }
 536   else{
 537     fCurrCentrality = -1;
 538   }
 539 }
 540
 541 //______________________________________________________________________________
 542 void AliPIDResponse::ExecNewRun()
 543 {
 544   //
 545   // Things to Execute upon a new run
 546   //
 547   SetRecoInfo();
 548
 549   SetITSParametrisation();
 550
 551   SetTPCPidResponseMaster();
 552   SetTPCParametrisation();
 553
 554   SetTRDPidResponseMaster();
 555   InitializeTRDResponse();
 556
 557   SetEMCALPidResponseMaster();
 558   InitializeEMCALResponse();
 559
 560   fTOFResponse.SetTimeResolution(fTOFres);
 561 }
 562
 563 //_____________________________________________________
 564 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
 565 {
 566   //
 567   // Get TPC multiplicity in bins of 150
 568   //
 569
 570   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
 571   Double_t tpcMulti=0.;
 572   if(vertexTPC){
 573     Double_t vertexContribTPC=vertexTPC->GetNContributors();
 574     tpcMulti=vertexContribTPC/150.;
 575     if (tpcMulti>20.) tpcMulti=20.;
 576   }
 577
 578   return tpcMulti;
 579 }
 580
 581 //______________________________________________________________________________
 582 void AliPIDResponse::SetRecoInfo()
 583 {
 584   //
 585   // Set reconstruction information
 586   //
 587
 588   //reset information
 589   fLHCperiod="";
 590   fMCperiodTPC="";
 591
 592   fBeamType="";
 593
 594   fBeamType="PP";
 595
 596   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
 597   //find the period by run number (UGLY, but not stored in ESD and AOD... )
 598   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
 599   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
 600   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
 601   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 602   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 603   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 604   else if (fRun>=136851&&fRun<=139517) {
 605     fLHCperiod="LHC10H";
 606     fMCperiodTPC="LHC10H8";
 607     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 608     fBeamType="PBPB";
 609   }
 610   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 611   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
 612   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 613   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 614   else if (fRun>=166529) {
 615     fLHCperiod="LHC11H";
 616     fMCperiodTPC="LHC11A10";
 617     fBeamType="PBPB";
 618   }
 619
 620
 621   //exception new pp MC productions from 2011
 622   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 623 }
 624
 625 //______________________________________________________________________________
 626 void AliPIDResponse::SetITSParametrisation()
 627 {
 628   //
 629   // Set the ITS parametrisation
 630   //
 631 }
 632
 633 //______________________________________________________________________________
 634 void AliPIDResponse::SetTPCPidResponseMaster()
 635 {
 636   //
 637   // Load the TPC pid response functions from the OADB
 638   //
 639   //don't load twice for the moment
 640    if (fArrPidResponseMaster) return;
 641
 642
 643   //reset the PID response functions
 644   delete fArrPidResponseMaster;
 645   fArrPidResponseMaster=0x0;
 646
 647   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
 648
 649   TFile *f=TFile::Open(fileName.Data());
 650   if (f && f->IsOpen() && !f->IsZombie()){
 651     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
 652   }
 653   delete f;
 654
 655   if (!fArrPidResponseMaster){
 656     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
 657     return;
 658   }
 659   fArrPidResponseMaster->SetOwner();
 660 }
 661
 662 //______________________________________________________________________________
 663 void AliPIDResponse::SetTPCParametrisation()
 664 {
 665   //
 666   // Change BB parametrisation for current run
 667   //
 668
 669   if (fLHCperiod.IsNull()) {
 670     AliFatal("No period set, not changing parametrisation");
 671     return;
 672   }
 673
 674   //
 675   // Set default parametrisations for data and MC
 676   //
 677
 678   //data type
 679   TString datatype="DATA";
 680   //in case of mc fRecoPass is per default 1
 681   if (fIsMC) {
 682     datatype="MC";
 683     fRecoPass=1;
 684   }
 685
 686   //
 687   //reset old splines
 688   //
 689   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 690     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
 691   }
 692
 693   //
 694   //set the new PID splines
 695   //
 696   TString period=fLHCperiod;
 697   if (fArrPidResponseMaster){
 698     TObject *grAll=0x0;
 699     //for MC don't use period information
 700 //     if (fIsMC) period="[A-Z0-9]*";
 701     //for MC use MC period information
 702     if (fIsMC) period=fMCperiodTPC;
 703 //pattern for the default entry (valid for all particles)
 704     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 705
 706     //loop over entries and filter them
 707     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
 708       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
 709       if (responseFunction==0x0) continue;
 710       TString responseName=responseFunction->GetName();
 711
 712       if (!reg.MatchB(responseName)) continue;
 713
 714       TObjArray *arr=reg.MatchS(responseName);
 715       TString particleName=arr->At(1)->GetName();
 716       delete arr;
 717       if (particleName.IsNull()) continue;
 718       if (particleName=="ALL") grAll=responseFunction;
 719       else {
 720         //find particle id
 721         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 722           TString particle=AliPID::ParticleName(ispec);
 723           particle.ToUpper();
 724           if ( particle == particleName ){
 725             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
 726             fTPCResponse.SetUseDatabase(kTRUE);
 727             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
 728             break;
 729           }
 730         }
 731       }
 732     }
 733
 734     //set default response function to all particles which don't have a specific one
 735     if (grAll){
 736       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 737         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
 738           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
 739           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
 740         }
 741       }
 742     }
 743   }
 744
 745   //
 746   // Setup resolution parametrisation
 747   //
 748
 749   //default
 750   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
 751
 752   if (fRun>=122195){
 753     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
 754   }
 755   if (fArrPidResponseMaster)
 756   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 757
 758   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
 759 }
 760
 761 //______________________________________________________________________________
 762 void AliPIDResponse::SetTRDPidResponseMaster()
 763 {
 764   //
 765   // Load the TRD pid params and references from the OADB
 766   //
 767   if(fTRDPIDParams) return;
 768   AliOADBContainer contParams("contParams");
 769
 770   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
 771   if(statusPars){
 772     AliError("Failed initializing PID Params from OADB");
 773   } else {
 774     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
 775     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
 776     if(!fTRDPIDParams){
 777       AliError(Form("TRD Params not found in run %d", fRun));
 778     }
 779   }
 780
 781   AliOADBContainer contRefs("contRefs");
 782   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
 783   if(statusRefs){
 784     AliInfo("Failed Loading References for TRD");
 785   } else {
 786     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
 787     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
 788     if(!fTRDPIDReference){
 789       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
 790     }
 791   }
 792 }
 793
 794 //______________________________________________________________________________
 795 void AliPIDResponse::InitializeTRDResponse(){
 796   //
 797   // Set PID Params and references to the TRD PID response
 798   //
 799   fTRDResponse.SetPIDParams(fTRDPIDParams);
 800   fTRDResponse.Load(fTRDPIDReference);
 801   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
 802     fTRDslicesForPID[0] = 0;
 803     fTRDslicesForPID[1] = 7;
 804   }
 805 }
 806
 807 //_________________________________________________________________________
 808 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
 809   //
 810   // Check whether track is identified as electron under a given electron efficiency hypothesis
 811   //
 812   Double_t probs[AliPID::kSPECIES];
 813   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
 814
 815   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
 816   // Take mean of the TRD momenta in the given tracklets
 817   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
 818   Int_t nmomenta = 0;
 819   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
 820     if(vtrack->GetTRDmomentum(iPl) > 0.){
 821       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
 822     }
 823   }
 824   p = TMath::Mean(nmomenta, trdmomenta);
 825
 826   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
 827 }
 828
 829 //______________________________________________________________________________
 830 void AliPIDResponse::SetEMCALPidResponseMaster()
 831 {
 832   //
 833   // Load the EMCAL pid response functions from the OADB
 834   //
 835   TObjArray* fEMCALPIDParamsRun      = NULL;
 836   TObjArray* fEMCALPIDParamsPass     = NULL;
 837
 838   if(fEMCALPIDParams) return;
 839   AliOADBContainer contParams("contParams");
 840
 841   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
 842   if(statusPars){
 843     AliError("Failed initializing PID Params from OADB");
 844   }
 845   else {
 846     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 847
 848     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
 849     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
 850     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 851
 852     if(!fEMCALPIDParams){
 853       AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
 854       AliInfo("Will take the standard LHC11d instead ...");
 855
 856       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
 857       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
 858       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 859
 860       if(!fEMCALPIDParams){
 861         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 862       }
 863     }
 864   }
 865 }
 866
 867 //______________________________________________________________________________
 868 void AliPIDResponse::InitializeEMCALResponse(){
 869   //
 870   // Set PID Params to the EMCAL PID response
 871   //
 872   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
 873
 874 }