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