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