PWGDQ/dielectron/AliDielectronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 ///////////////////////////////////////////////////////////////////////////
  17 //                Dielectron PID                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Detailed description
  22
  23
  24 */
  25 //                                                                       //
  26 ///////////////////////////////////////////////////////////////////////////
  27
  28 #include <TMath.h>
  29 #include <TF1.h>
  30 #include <TGraph.h>
  31 #include <THnBase.h>
  32
  33 #include <AliVTrack.h>
  34 #include <AliVCluster.h>
  35 #include <AliLog.h>
  36 #include <AliExternalTrackParam.h>
  37 #include <AliPIDResponse.h>
  38 #include <AliTRDPIDResponse.h>
  39 #include <AliESDtrack.h> //!!!!! Remove once Eta correction is treated in the tender
  40 #include <AliAODTrack.h>
  41 #include <AliAODPid.h>
  42
  43 #include "AliDielectronVarManager.h"
  44 #include "AliDielectronVarCuts.h"
  45
  46 #include "AliDielectronPID.h"
  47
  48 ClassImp(AliDielectronPID)
  49
  50 TGraph  *AliDielectronPID::fgFitCorr=0x0;
  51 Double_t AliDielectronPID::fgCorr=0.0;
  52 Double_t AliDielectronPID::fgCorrdEdx=1.0;
  53 TF1     *AliDielectronPID::fgFunEtaCorr=0x0;
  54 TH1     *AliDielectronPID::fgFunCntrdCorr=0x0;
  55 TH1     *AliDielectronPID::fgFunWdthCorr=0x0;
  56 TGraph  *AliDielectronPID::fgdEdxRunCorr=0x0;
  57
  58 AliDielectronPID::AliDielectronPID() :
  59   AliAnalysisCuts(),
  60   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  61   fNcuts(0),
  62   fPIDResponse(0x0)
  63 {
  64   //
  65   // Default Constructor
  66   //
  67   for (Int_t icut=0; icut<kNmaxPID; ++icut){
  68     fDetType[icut]=kTPC;
  69     fPartType[icut]=AliPID::kPion;
  70     fNsigmaLow[icut]=0;
  71     fNsigmaUp[icut]=0;
  72     fmin[icut]=0;
  73     fmax[icut]=0;
  74     fExclude[icut]=kFALSE;
  75     fFunUpperCut[icut]=0x0;
  76     fFunLowerCut[icut]=0x0;
  77     fRequirePIDbit[icut]=0;
  78     fActiveCuts[icut]=-1;
  79     fSigmaFunLow[icut]=0;
  80     fSigmaFunUp[icut]=0;
  81     fFunSigma[icut]=0x0;
  82     fVarCuts[icut]=0x0;
  83     fMapElectronCutLow[icut]=0x0;
  84   }
  85 }
  86
  87 //______________________________________________
  88 AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
  89   AliAnalysisCuts(name, title),
  90   fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
  91   fNcuts(0),
  92   fPIDResponse(0x0)
  93 {
  94   //
  95   // Named Constructor
  96   //
  97   for (Int_t icut=0; icut<kNmaxPID; ++icut){
  98     fDetType[icut]=kTPC;
  99     fPartType[icut]=AliPID::kPion;
 100     fNsigmaLow[icut]=0;
 101     fNsigmaUp[icut]=0;
 102     fmin[icut]=0;
 103     fmax[icut]=0;
 104     fExclude[icut]=kFALSE;
 105     fFunUpperCut[icut]=0x0;
 106     fFunLowerCut[icut]=0x0;
 107     fRequirePIDbit[icut]=0;
 108     fActiveCuts[icut]=0;
 109     fSigmaFunLow[icut]=0;
 110     fSigmaFunUp[icut]=0;
 111     fFunSigma[icut]=0x0;
 112     fVarCuts[icut]=0x0;
 113     fMapElectronCutLow[icut]=0x0;
 114   }
 115 }
 116
 117 //______________________________________________
 118 AliDielectronPID::~AliDielectronPID()
 119 {
 120   //
 121   // Default Destructor
 122   //
 123   if (fUsedVars) delete fUsedVars;
 124 }
 125
 126 //______________________________________________
 127 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/*=-99999.*/,
 128                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 129                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 130 {
 131   //
 132   // Add a pid nsigma cut
 133   // use response of detector 'det' in the range ['min'] to ['max'] for var
 134   // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
 135   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 136   // specify whether to 'exclude' the given band
 137   //
 138
 139   if (fNcuts>=kNmaxPID){
 140     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 141     return;
 142   }
 143   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 144     nSigmaUp=TMath::Abs(nSigmaLow);
 145     nSigmaLow=-1.*nSigmaUp;
 146   }
 147   fDetType[fNcuts]=det;
 148   fPartType[fNcuts]=type;
 149   fNsigmaLow[fNcuts]=nSigmaLow;
 150   fNsigmaUp[fNcuts]=nSigmaUp;
 151   fmin[fNcuts]=min;
 152   fmax[fNcuts]=max;
 153   fExclude[fNcuts]=exclude;
 154   fRequirePIDbit[fNcuts]=pidBitType;
 155   if(var==-1) var=AliDielectronVarManager::kP; // set default
 156   fActiveCuts[fNcuts]=var;
 157   fUsedVars->SetBitNumber(var,kTRUE);
 158
 159   AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
 160                   fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));
 161
 162   ++fNcuts;
 163
 164 }
 165
 166 //______________________________________________
 167 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
 168                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 169                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 170 {
 171   //
 172   // cut using a TF1 as upper cut
 173   //
 174   if (funUp==0x0){
 175     AliError("A valid function is required for the upper cut. Not adding the cut!");
 176     return;
 177   }
 178   fFunUpperCut[fNcuts]=funUp;
 179   AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var);
 180 }
 181
 182 //______________________________________________
 183 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
 184                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 185                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 186 {
 187   //
 188   // cut using a TF1 as lower cut
 189   //
 190   if (funLow==0x0){
 191     AliError("A valid function is required for the lower cut. Not adding the cut!");
 192     return;
 193   }
 194   fFunLowerCut[fNcuts]=funLow;
 195   AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
 196 }
 197
 198 //______________________________________________
 199 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
 200                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 201                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 202 {
 203   //
 204   // cut using a TF1 as lower and upper cut
 205   //
 206   if ( (funUp==0x0) || (funLow==0x0) ){
 207     AliError("A valid function is required for upper and lower cut. Not adding the cut!");
 208     return;
 209   }
 210   fFunUpperCut[fNcuts]=funUp;
 211   fFunLowerCut[fNcuts]=funLow;
 212   AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
 213 }
 214
 215 //______________________________________________
 216 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
 217                               Double_t min, Double_t max, Bool_t exclude,
 218                               UInt_t pidBitType, TF1 * const funSigma)
 219 {
 220   //
 221   // cut using a TF1 as lower cut
 222   //
 223   if (funSigma==0x0){
 224     AliError("A valid function is required for the sigma cut. Not adding the cut!");
 225     return;
 226   }
 227   fFunSigma[fNcuts]=funSigma;
 228   fSigmaFunLow[fNcuts]=min;
 229   fSigmaFunUp[fNcuts]=max;
 230
 231   AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1);
 232 }
 233
 234 //______________________________________________
 235 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, THnBase * const histLow, Double_t nSigmaUp,
 236                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 237                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 238 {
 239   //
 240   // cut using a THnBase as a lower cut
 241   //
 242   if(histLow==0x0){
 243     AliError("A valid histogram is required for the lower cut. Not adding the cut!");
 244     return;
 245   }
 246   fMapElectronCutLow[fNcuts]=histLow;
 247   AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
 248 }
 249 //______________________________________________
 250 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
 251                               AliDielectronVarCuts *var, Bool_t exclude/*=kFALSE*/,
 252                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 253 {
 254   //
 255   // Add a pid nsigma cut
 256   // use response of detector 'det' in the ranges for variables defined in var
 257   // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
 258   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 259   // specify whether to 'exclude' the given band
 260   //
 261   if(!var) return;
 262   if (fNcuts>=kNmaxPID){
 263     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 264     return;
 265   }
 266   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 267     nSigmaUp=TMath::Abs(nSigmaLow);
 268     nSigmaLow=-1.*nSigmaUp;
 269   }
 270   fDetType[fNcuts]=det;
 271   fPartType[fNcuts]=type;
 272   fNsigmaLow[fNcuts]=nSigmaLow;
 273   fNsigmaUp[fNcuts]=nSigmaUp;
 274   fExclude[fNcuts]=exclude;
 275   fRequirePIDbit[fNcuts]=pidBitType;
 276   fVarCuts[fNcuts]=var;
 277
 278   AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \n",
 279                   fNcuts,nSigmaLow,nSigmaUp));
 280
 281   ++fNcuts;
 282
 283 }
 284
 285 //______________________________________________
 286 Bool_t AliDielectronPID::IsSelected(TObject* track)
 287 {
 288   //
 289   // perform PID cuts
 290   //
 291
 292   //loop over all cuts
 293   AliVTrack *part=static_cast<AliVTrack*>(track);
 294   AliESDtrack *esdTrack=0x0;
 295   AliAODTrack *aodTrack=0x0;
 296   Double_t origdEdx=-1;
 297
 298   // apply ETa correction, remove once this is in the tender
 299   if( (part->IsA() == AliESDtrack::Class()) ){
 300     esdTrack=static_cast<AliESDtrack*>(part);
 301     origdEdx=esdTrack->GetTPCsignal();
 302     esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 303   } else if ( (part->IsA() == AliAODTrack::Class()) ){
 304     aodTrack=static_cast<AliAODTrack*>(track);
 305     AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 306     if (pid){
 307       origdEdx=pid->GetTPCsignal();
 308       pid->SetTPCsignal(origdEdx/GetEtaCorr(aodTrack)/fgCorrdEdx);
 309     }
 310   }
 311
 312   // check for corrections and add their variables to the fill map
 313   if(fgFunCntrdCorr)  {
 314     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetXaxis()->GetUniqueID(), kTRUE);
 315     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetYaxis()->GetUniqueID(), kTRUE);
 316     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetZaxis()->GetUniqueID(), kTRUE);
 317   }
 318   if(fgFunWdthCorr)  {
 319     fUsedVars->SetBitNumber(fgFunWdthCorr->GetXaxis()->GetUniqueID(), kTRUE);
 320     fUsedVars->SetBitNumber(fgFunWdthCorr->GetYaxis()->GetUniqueID(), kTRUE);
 321     fUsedVars->SetBitNumber(fgFunWdthCorr->GetZaxis()->GetUniqueID(), kTRUE);
 322   }
 323   for(UChar_t icut=0; icut<fNcuts; ++icut) {
 324     if(!fMapElectronCutLow[icut]) continue;
 325     for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
 326       TString var(fMapElectronCutLow[icut]->GetAxis(idim)->GetName());
 327       fUsedVars->SetBitNumber(AliDielectronVarManager::GetValueType(var.Data()), kTRUE);
 328     }
 329   }
 330
 331   //Fill values
 332   Double_t values[AliDielectronVarManager::kNMaxValues];
 333   AliDielectronVarManager::SetFillMap(fUsedVars);
 334   AliDielectronVarManager::Fill(track,values);
 335
 336   Bool_t selected=kFALSE;
 337   fPIDResponse=AliDielectronVarManager::GetPIDResponse();
 338   for (UChar_t icut=0; icut<fNcuts; ++icut){
 339     Double_t min=fmin[icut];
 340     Double_t max=fmax[icut];
 341     Double_t val=values[fActiveCuts[icut]];
 342
 343     // test var range. In case min==max do not cut
 344     if ( fVarCuts[icut] ) {
 345       if ( !fVarCuts[icut]->IsSelected(part) ) continue;
 346     }
 347     else if ( ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) ) {
 348         continue;
 349     }
 350
 351     // check if fFunSigma is set, then check if 'part' is in sigma range of the function
 352     if(fFunSigma[icut]){
 353         val= fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
 354         if (fPartType[icut]==AliPID::kElectron){
 355             val-=fgCorr;
 356         }
 357         min= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunLow[icut];
 358         max= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunUp[icut];
 359         if(val<min || val>=max) continue;
 360     }
 361
 362     switch (fDetType[icut]){
 363     case kITS:
 364       selected = IsSelectedITS(part,icut);
 365       break;
 366     case kTPC:
 367       selected = IsSelectedTPC(part,icut,values);
 368       break;
 369     case kTRD:
 370       selected = IsSelectedTRD(part,icut);
 371       break;
 372     case kTRDeleEff:
 373       selected = IsSelectedTRDeleEff(part,icut);
 374       break;
 375     case kTRDeleEff2D:
 376       selected = IsSelectedTRDeleEff(part,icut,AliTRDPIDResponse::kLQ2D);
 377       break;
 378     case kTOF:
 379       selected = IsSelectedTOF(part,icut);
 380       break;
 381     case kEMCAL:
 382       selected = IsSelectedEMCAL(part,icut);
 383       break;
 384     }
 385     if (!selected) {
 386       if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 387       else if (aodTrack){
 388         AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 389         if (pid) pid->SetTPCsignal(origdEdx);
 390       }
 391
 392       return kFALSE;
 393     }
 394   }
 395
 396   if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 397   else if (aodTrack){
 398     AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 399     if (pid) pid->SetTPCsignal(origdEdx);
 400   }
 401   return selected;
 402 }
 403
 404 //______________________________________________
 405 Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
 406 {
 407   //
 408   // ITS part of the PID check
 409   // Don't accept the track if there was no pid bit set
 410   //
 411   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kITS,part);
 412   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 413   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 414
 415   Double_t mom=part->P();
 416
 417   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]);
 418
 419   // test if we are supposed to use a function for the cut
 420   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
 421   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
 422
 423   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 424   return selected;
 425 }
 426
 427 //______________________________________________
 428 Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut, Double_t *values)
 429 {
 430   //
 431   // TPC part of the PID check
 432   // Don't accept the track if there was no pid bit set
 433   //
 434   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,part);
 435   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 436   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 437
 438
 439   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
 440
 441   // eta corrections
 442   if (fPartType[icut]==AliPID::kElectron){
 443     // old way 1D
 444     numberOfSigmas-=fgCorr;
 445     // via functions (1-3D)
 446     numberOfSigmas-=GetCntrdCorr(part);
 447     numberOfSigmas/=GetWdthCorr(part);
 448   }
 449
 450   // matching of MC and data //
 451
 452   // test if we are supposed to use a function for the cut (matching via fits)
 453   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
 454   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
 455
 456   // test if we are using cut maps (in calibrated units)
 457   Double_t lowElectronCut = fNsigmaLow[icut];
 458   Double_t upElectronCut  = fNsigmaUp[icut];
 459
 460   if((fPartType[icut]==AliPID::kElectron) && (fMapElectronCutLow[icut] )) {
 461     Double_t *vals = new Double_t[fMapElectronCutLow[icut]->GetNdimensions()];//={-1};
 462     // get array of values for the corresponding dimensions using axis names
 463     for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
 464       vals[idim] = values[AliDielectronVarManager::GetValueType(fMapElectronCutLow[icut]->GetAxis(idim)->GetName())];
 465       // printf(" \t %s %.3f ",fMapElectronCutLow[icut]->GetAxis(idim)->GetName(),vals[idim]);
 466     }
 467     // find bin for values (w/o creating it in case it is not filled)
 468     Long_t bin = fMapElectronCutLow[icut]->GetBin(vals,kFALSE);
 469     if(bin>0) lowElectronCut=fMapElectronCutLow[icut]->GetBinContent(bin);
 470     else lowElectronCut=100;
 471     //    printf("  low cut %.3f (%ld) \n", lowElectronCut,bin);
 472     delete [] vals;
 473   }
 474
 475
 476   Bool_t selected=((numberOfSigmas>=lowElectronCut)&&(numberOfSigmas<=upElectronCut))^fExclude[icut];
 477   return selected;
 478 }
 479
 480 //______________________________________________
 481 Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
 482 {
 483   //
 484   // TRD part of the pid check
 485   // the TRD checks on the probabilities.
 486   //
 487
 488   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 489   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 490   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 491
 492   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
 493
 494   Double_t p[AliPID::kSPECIES]={0.};
 495   fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
 496   Float_t particleProb=p[fPartType[icut]];
 497
 498   Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
 499   return selected;
 500 }
 501
 502 //______________________________________________
 503 Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod)
 504 {
 505   //
 506   // TRD part of the pid check using electron efficiency requirement
 507   // in this case the upper limit as well as the particle specie is ignored
 508   //   and the lower limit regarded as the requested electron efficiency
 509   //
 510
 511   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 512   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 513   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 514
 515   Double_t centrality = -1.;
 516   if(part->IsA() == AliESDtrack::Class())
 517     centrality=(const_cast<AliESDEvent*>( (static_cast<const AliESDtrack*>(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 518   if(part->IsA() == AliAODTrack::Class())
 519     centrality=(const_cast<AliAODEvent*>( (static_cast<const AliAODTrack*>(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 520
 521   Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut];
 522   return selected;
 523 }
 524
 525 //______________________________________________
 526 Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
 527 {
 528   //
 529   // TOF part of the PID check
 530   // Don't accept the track if there was no pid bit set
 531   //
 532   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,part);
 533   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 534   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 535
 536   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
 537
 538   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 539   return selected;
 540 }
 541
 542 //______________________________________________
 543 Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
 544 {
 545   //
 546   // emcal pid selecttion
 547   //
 548
 549   //TODO: correct way to check for emcal pid?
 550   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);
 551
 552   Bool_t hasPID=numberOfSigmas>-998.;
 553
 554   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
 555   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;
 556
 557
 558   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 559   return selected;
 560 }
 561
 562 //______________________________________________
 563 void AliDielectronPID::SetDefaults(Int_t def){
 564   //
 565   // initialise default pid strategies
 566   //
 567
 568   if (def==0){
 569     // 2sigma bands TPC:
 570     // - include e
 571     // - exclude mu,K,pi,p
 572     // -complete p range
 573     AddCut(kTPC,AliPID::kElectron,2);
 574     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
 575     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
 576     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
 577     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
 578   } else if (def==1) {
 579     // 2sigma bands TPC:
 580     // - include e 0<p<inf
 581     // - exclude mu,K,pi,p 0<p<2
 582     AddCut(kTPC,AliPID::kElectron,2.);
 583     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
 584     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
 585     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
 586     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
 587
 588   } else if (def==2) {
 589     // include 2sigma e TPC
 590     // 3sigma bands TOF
 591     // - exclude K,P
 592     AddCut(kTPC,AliPID::kElectron,-3.,3.);
 593     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 594     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 595
 596   } else if (def==3 || def==4) { // def3 and def4 are the same??
 597     // include 2sigma e TPC
 598     // 3sigma bands TOF
 599     // - exclude K 0<p<1
 600     // - exclude P 0<p<2
 601     AddCut(kTPC,AliPID::kElectron,2);
 602     AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
 603     AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
 604     AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
 605
 606   } else if (def==5) {
 607     AddCut(kTPC,AliPID::kElectron,-0.5,3);
 608     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 609   } else if (def==6) {
 610     // lower cut TPC: parametrisation by HFE
 611     // upper cut TPC: 3 sigma
 612     // TOF ele band 3sigma 0<p<1.5GeV
 613     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 614     lowerCut->SetParameters(-2.7,-0.4357);
 615     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 616     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 617   } else if (def==7) {
 618     // wide TPC cut
 619     // TOF ele band 3sigma 0<p<1.5GeV
 620     AddCut(kTPC,AliPID::kElectron,10.);
 621     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 622   } else if (def==8) {
 623     // TOF 5 sigma inclusion if TOFpid available
 624     // this should reduce K,p,Pi to a large extent
 625     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 626   } else if (def==9) {
 627     // lower cut TPC: parametrisation by HFE
 628     // upper cut TPC: 3 sigma
 629     // TOF 5 sigma inclusion if TOFpid available
 630     // this should reduce K,p,Pi to a large extent
 631     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 632     lowerCut->SetParameters(-2.65,-0.6757);
 633     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 634     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 635
 636   } else if (def==10) {
 637     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 638     AddCut(kTPC,AliPID::kElectron,3.);
 639     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 640     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 641
 642   } else if (def==11) {
 643     // lower cut TPC: parametrisation by HFE
 644     // only use from period d on !!
 645     // upper cut TPC: 3 sigma
 646     // TOF ele band 3sigma 0<p<2.0GeV
 647     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 648     lowerCut->SetParameters(-3.718,-0.4,0.27);
 649     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 650     AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
 651   } else if (def==12) {
 652     // lower cut TPC: parametrisation by HFE
 653     // only use from period d on !!
 654     // upper cut TPC: 3 sigma
 655     // TOF 5 sigma inclusion if TOFpid available
 656     // this should reduce K,p,Pi to a large extent
 657     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 658     lowerCut->SetParameters(-3.718,-0.4,0.27);
 659     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 660     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 661   }
 662   else if (def==13) {
 663     // TPC electron inclusion
 664     // TOF electron inclusion if available
 665     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 666     AddCut(kTPC,AliPID::kElectron,10.);
 667   }
 668   else if (def==14) {
 669     // TRD 1D 90% elec eff, 4-6 tracklets
 670     // TPC electron inclusion
 671     // TOF electron inclusion if available
 672     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 673            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 674     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 675     AddCut(kTPC,AliPID::kElectron,10.);
 676   }
 677   else if (def==15) {
 678     // TRD 1D 90% elec eff, 4-6 tracklets, chi2 < 2
 679     // TPC electron inclusion
 680     // TOF electron inclusion if available
 681     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 682            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 683     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,0.,2.,kFALSE,
 684            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDchi2);
 685     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 686     AddCut(kTPC,AliPID::kElectron,10.);
 687   }
 688   else if (def==16) {
 689     // TPC electron inclusion
 690     // TOF electron inclusion if available
 691     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 692     AddCut(kTPC,AliPID::kElectron,-3.5,+3.5);
 693   }
 694
 695 }
 696
 697 //______________________________________________
 698 void AliDielectronPID::SetCorrVal(Double_t run)
 699 {
 700   //
 701   // set correction value for run
 702   //
 703   fgCorr=0.;
 704   fgCorrdEdx=1.;
 705
 706   if (fgFitCorr){
 707     fgCorr=fgFitCorr->Eval(run);
 708     if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
 709     if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
 710   }
 711
 712   if (fgdEdxRunCorr){
 713     fgCorrdEdx=fgdEdxRunCorr->Eval(run);
 714     if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
 715     if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
 716   }
 717 }
 718
 719 //______________________________________________
 720 Double_t AliDielectronPID::GetEtaCorr(const AliVTrack *track)
 721 {
 722   //
 723   // return eta correction
 724   //
 725   if (!fgFunEtaCorr) return 1;
 726   return fgFunEtaCorr->Eval(track->Eta());
 727 }
 728
 729 //______________________________________________
 730 Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TH1 *hist)
 731 {
 732   //
 733   // return correction value
 734   //
 735   //TODO: think about taking an values array as argument to reduce # var fills
 736
 737   //Fill only event and vparticle values (otherwise we end up in a circle)
 738   Double_t values[AliDielectronVarManager::kNMaxValues];
 739   AliDielectronVarManager::FillVarVParticle(track,values);
 740
 741   TF1 *fun = (TF1*)hist->GetListOfFunctions()->At(0);
 742   Int_t dim=fun->GetNdim();
 743
 744   Double_t var[3] = {0.,0.,0.};
 745   if(dim>0) var[0] = values[hist->GetXaxis()->GetUniqueID()];
 746   if(dim>1) var[1] = values[hist->GetYaxis()->GetUniqueID()];
 747   if(dim>2) var[2] = values[hist->GetZaxis()->GetUniqueID()];
 748   Double_t corr = fun->Eval(var[0],var[1],var[2]);
 749   //  printf("%d-dim CORR value: %f (track %p) \n",dim,corr,track);
 750   return corr;
 751 }