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