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 <TH2D.h>
  32 #include <TH3D.h>
  33
  34 #include <AliVTrack.h>
  35 #include <AliVCluster.h>
  36 #include <AliLog.h>
  37 #include <AliExternalTrackParam.h>
  38 #include <AliPIDResponse.h>
  39 #include <AliTRDPIDResponse.h>
  40 #include <AliESDtrack.h> //!!!!! Remove once Eta correction is treated in the tender
  41 #include <AliAODTrack.h>
  42 #include <AliAODPid.h>
  43
  44 #include "AliDielectronVarManager.h"
  45 #include "AliDielectronVarCuts.h"
  46
  47 #include "AliDielectronPID.h"
  48
  49 ClassImp(AliDielectronPID)
  50
  51 TGraph *AliDielectronPID::fgFitCorr=0x0;
  52 Double_t AliDielectronPID::fgCorr=0.0;
  53 Double_t AliDielectronPID::fgCorrdEdx=1.0;
  54 TF1 *AliDielectronPID::fgFunEtaCorr=0x0;
  55 TF1 *AliDielectronPID::fgFunCntrdCorr=0x0;
  56 TF1 *AliDielectronPID::fgFunWdthCorr=0x0;
  57 TGraph *AliDielectronPID::fgdEdxRunCorr=0x0;
  58
  59 AliDielectronPID::AliDielectronPID() :
  60   AliAnalysisCuts(),
  61   fNcuts(0),
  62   fPIDResponse(0x0),
  63   fElectronCentroidCentEta(0x0),
  64   fElectronWidthCentEta(0x0)
  65 {
  66   //
  67   // Default Constructor
  68   //
  69   for (Int_t icut=0; icut<kNmaxPID; ++icut){
  70     fDetType[icut]=kTPC;
  71     fPartType[icut]=AliPID::kPion;
  72     fNsigmaLow[icut]=0;
  73     fNsigmaUp[icut]=0;
  74     fmin[icut]=0;
  75     fmax[icut]=0;
  76     fExclude[icut]=kFALSE;
  77     fFunUpperCut[icut]=0x0;
  78     fFunLowerCut[icut]=0x0;
  79     fRequirePIDbit[icut]=0;
  80     fActiveCuts[icut]=-1;
  81     fSigmaFunLow[icut]=0;
  82     fSigmaFunUp[icut]=0;
  83     fFunSigma[icut]=0x0;
  84     fVarCuts[icut]=0x0;
  85     fHistElectronCutLow[icut]=0x0;
  86     fHistElectronCutUp[icut]=0x0;
  87   }
  88 }
  89
  90 //______________________________________________
  91 AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
  92   AliAnalysisCuts(name, title),
  93   fNcuts(0),
  94   fPIDResponse(0x0),
  95   fElectronCentroidCentEta(0x0),
  96   fElectronWidthCentEta(0x0)
  97 {
  98   //
  99   // Named Constructor
 100   //
 101   for (Int_t icut=0; icut<kNmaxPID; ++icut){
 102     fDetType[icut]=kTPC;
 103     fPartType[icut]=AliPID::kPion;
 104     fNsigmaLow[icut]=0;
 105     fNsigmaUp[icut]=0;
 106     fmin[icut]=0;
 107     fmax[icut]=0;
 108     fExclude[icut]=kFALSE;
 109     fFunUpperCut[icut]=0x0;
 110     fFunLowerCut[icut]=0x0;
 111     fRequirePIDbit[icut]=0;
 112     fActiveCuts[icut]=0;
 113     fSigmaFunLow[icut]=0;
 114     fSigmaFunUp[icut]=0;
 115     fFunSigma[icut]=0x0;
 116     fVarCuts[icut]=0x0;
 117     fHistElectronCutLow[icut]=0x0;
 118     fHistElectronCutUp[icut]=0x0;
 119   }
 120 }
 121
 122 //______________________________________________
 123 AliDielectronPID::~AliDielectronPID()
 124 {
 125   //
 126   // Default Destructor
 127   //
 128 }
 129
 130 //______________________________________________
 131 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/*=-99999.*/,
 132                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 133                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 134 {
 135   //
 136   // Add a pid nsigma cut
 137   // use response of detector 'det' in the range ['min'] to ['max'] for var
 138   // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
 139   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 140   // specify whether to 'exclude' the given band
 141   //
 142
 143   if (fNcuts>=kNmaxPID){
 144     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 145     return;
 146   }
 147   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 148     nSigmaUp=TMath::Abs(nSigmaLow);
 149     nSigmaLow=-1.*nSigmaUp;
 150   }
 151   fDetType[fNcuts]=det;
 152   fPartType[fNcuts]=type;
 153   fNsigmaLow[fNcuts]=nSigmaLow;
 154   fNsigmaUp[fNcuts]=nSigmaUp;
 155   fmin[fNcuts]=min;
 156   fmax[fNcuts]=max;
 157   fExclude[fNcuts]=exclude;
 158   fRequirePIDbit[fNcuts]=pidBitType;
 159   fActiveCuts[fNcuts]=(var==-1 ? AliDielectronVarManager::kP : var);
 160
 161   AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n",
 162                   fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts])));
 163
 164   ++fNcuts;
 165
 166 }
 167
 168 //______________________________________________
 169 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
 170                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 171                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 172 {
 173   //
 174   // cut using a TF1 as upper cut
 175   //
 176   if (funUp==0x0){
 177     AliError("A valid function is required for the upper cut. Not adding the cut!");
 178     return;
 179   }
 180   fFunUpperCut[fNcuts]=funUp;
 181   AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var);
 182 }
 183
 184 //______________________________________________
 185 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
 186                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 187                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 188 {
 189   //
 190   // cut using a TF1 as lower cut
 191   //
 192   if (funLow==0x0){
 193     AliError("A valid function is required for the lower cut. Not adding the cut!");
 194     return;
 195   }
 196   fFunLowerCut[fNcuts]=funLow;
 197   AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
 198 }
 199
 200 //______________________________________________
 201 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
 202                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 203                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 204 {
 205   //
 206   // cut using a TF1 as lower and upper cut
 207   //
 208   if ( (funUp==0x0) || (funLow==0x0) ){
 209     AliError("A valid function is required for upper and lower cut. Not adding the cut!");
 210     return;
 211   }
 212   fFunUpperCut[fNcuts]=funUp;
 213   fFunLowerCut[fNcuts]=funLow;
 214   AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
 215 }
 216
 217 //______________________________________________
 218 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
 219                               Double_t min, Double_t max, Bool_t exclude,
 220                               UInt_t pidBitType, TF1 * const funSigma)
 221 {
 222   //
 223   // cut using a TF1 as lower cut
 224   //
 225   if (funSigma==0x0){
 226     AliError("A valid function is required for the sigma cut. Not adding the cut!");
 227     return;
 228   }
 229   fFunSigma[fNcuts]=funSigma;
 230   fSigmaFunLow[fNcuts]=min;
 231   fSigmaFunUp[fNcuts]=max;
 232
 233   AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1);
 234 }
 235
 236 //______________________________________________
 237 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TH3D * const histLow, Double_t nSigmaUp,
 238                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 239                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 240 {
 241   //
 242   // cut using a TH3D(Pin,centrality,eta) as a lower cut
 243   //
 244   if(histLow==0x0){
 245     AliError("A valid histogram is required for the lower cut. Not adding the cut!");
 246     return;
 247   }
 248   fHistElectronCutLow[fNcuts]=histLow;
 249   AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
 250 }
 251
 252 //______________________________________________
 253 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TH3D * const histLow, TH3D * const histUp,
 254                               Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/,
 255                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/)
 256 {
 257   //
 258   // cut using a TH3D(Pin,centrality,eta) as lower and upper cut
 259   //
 260   if ( (histUp==0x0) || (histLow==0x0) ){
 261     AliError("A valid histogram is required for upper and lower cut. Not adding the cut!");
 262     return;
 263   }
 264   fHistElectronCutUp[fNcuts]=histUp;
 265   fHistElectronCutLow[fNcuts]=histLow;
 266   AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
 267 }
 268
 269 //______________________________________________
 270 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
 271                               AliDielectronVarCuts *var, Bool_t exclude/*=kFALSE*/,
 272                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 273 {
 274   //
 275   // Add a pid nsigma cut
 276   // use response of detector 'det' in the ranges for variables defined in var
 277   // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
 278   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 279   // specify whether to 'exclude' the given band
 280   //
 281   if(!var) return;
 282   if (fNcuts>=kNmaxPID){
 283     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 284     return;
 285   }
 286   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 287     nSigmaUp=TMath::Abs(nSigmaLow);
 288     nSigmaLow=-1.*nSigmaUp;
 289   }
 290   fDetType[fNcuts]=det;
 291   fPartType[fNcuts]=type;
 292   fNsigmaLow[fNcuts]=nSigmaLow;
 293   fNsigmaUp[fNcuts]=nSigmaUp;
 294   fExclude[fNcuts]=exclude;
 295   fRequirePIDbit[fNcuts]=pidBitType;
 296   fVarCuts[fNcuts]=var;
 297
 298   AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \n",
 299                   fNcuts,nSigmaLow,nSigmaUp));
 300
 301   ++fNcuts;
 302
 303 }
 304
 305 //______________________________________________
 306 Bool_t AliDielectronPID::IsSelected(TObject* track)
 307 {
 308   //
 309   // perform PID cuts
 310   //
 311
 312   //loop over all cuts
 313   AliVTrack *part=static_cast<AliVTrack*>(track);
 314   AliESDtrack *esdTrack=0x0;
 315   AliAODTrack *aodTrack=0x0;
 316   Double_t origdEdx=-1;
 317
 318   // apply ETa correction, remove once this is in the tender
 319   if( (part->IsA() == AliESDtrack::Class()) ){
 320     esdTrack=static_cast<AliESDtrack*>(part);
 321     origdEdx=esdTrack->GetTPCsignal();
 322     esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 323   } else if ( (part->IsA() == AliAODTrack::Class()) ){
 324     aodTrack=static_cast<AliAODTrack*>(track);
 325     AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 326     if (pid){
 327       origdEdx=pid->GetTPCsignal();
 328       pid->SetTPCsignal(origdEdx/GetEtaCorr(aodTrack)/fgCorrdEdx);
 329     }
 330   }
 331
 332   //Fill values
 333   Double_t values[AliDielectronVarManager::kNMaxValues];
 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);
 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)
 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   Double_t mom=part->GetTPCmomentum();
 439   Double_t centroid=0.0;
 440   Double_t width=1.0;
 441   Double_t eta = 0.0;
 442   Double_t centrality = 0.0;
 443
 444   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
 445
 446   if (fPartType[icut]==AliPID::kElectron){
 447     numberOfSigmas-=fgCorr;
 448     numberOfSigmas-=GetCntrdCorr(part);
 449     numberOfSigmas/=GetWdthCorr(part);
 450
 451     // recalculate the number of sigmas from the electron assuption using the electron centroid and width maps
 452     if(fElectronCentroidCentEta || fElectronWidthCentEta) {
 453       eta=part->Eta();
 454       centrality = AliDielectronVarManager::GetCurrentEvent()->GetCentrality()->GetCentralityPercentile("V0M");
 455       if(fElectronCentroidCentEta) {
 456         Int_t centBin = fElectronCentroidCentEta->GetXaxis()->FindBin(centrality);
 457         Int_t etaBin = fElectronCentroidCentEta->GetYaxis()->FindBin(eta);
 458         if(centBin>=1 && centBin<=fElectronCentroidCentEta->GetXaxis()->GetNbins() &&
 459            etaBin>=1 && etaBin<=fElectronCentroidCentEta->GetYaxis()->GetNbins())
 460           centroid = fElectronCentroidCentEta->GetBinContent(fElectronCentroidCentEta->GetXaxis()->FindBin(centrality),
 461                                                              fElectronCentroidCentEta->GetYaxis()->FindBin(eta));
 462       }
 463       if(fElectronWidthCentEta) {
 464         Int_t centBin =fElectronWidthCentEta->GetXaxis()->FindBin(centrality);
 465         Int_t etaBin = fElectronWidthCentEta->GetYaxis()->FindBin(eta);
 466         if(centBin>=1 && centBin<=fElectronWidthCentEta->GetXaxis()->GetNbins() &&
 467            etaBin>=1 && etaBin<=fElectronWidthCentEta->GetYaxis()->GetNbins())
 468           width = fElectronWidthCentEta->GetBinContent(fElectronWidthCentEta->GetXaxis()->FindBin(centrality),
 469                                                        fElectronWidthCentEta->GetYaxis()->FindBin(eta));
 470       }
 471     }
 472     if(width<0.01) width = 1.0;
 473     numberOfSigmas = (numberOfSigmas-centroid)/width;
 474   }
 475
 476   // test if we are supposed to use a function for the cut
 477   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
 478   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
 479
 480   // if electron selection 3D maps are loaded, then find the corresponding cuts based on momentum, eta and centrality
 481   Double_t lowElectronCut = fNsigmaLow[icut];
 482   Double_t upElectronCut = fNsigmaUp[icut];
 483   if((fPartType[icut]==AliPID::kElectron) && (fHistElectronCutLow[icut] || fHistElectronCutUp[icut])) {
 484     eta=part->Eta();
 485     centrality = AliDielectronVarManager::GetCurrentEvent()->GetCentrality()->GetCentralityPercentile("V0M");
 486   }
 487   if((fPartType[icut]==AliPID::kElectron) && fHistElectronCutLow[icut]) {
 488     Int_t pbin = fHistElectronCutLow[icut]->GetZaxis()->FindBin(mom);
 489     Int_t centBin = fHistElectronCutLow[icut]->GetXaxis()->FindBin(centrality);
 490     Int_t etaBin = fHistElectronCutLow[icut]->GetYaxis()->FindBin(eta);
 491     // apply this cut only in the range covered in the histogram range
 492     if(pbin!=0    && pbin!=(fHistElectronCutLow[icut]->GetZaxis()->GetNbins()+1)             &&
 493        centBin!=0 && centBin!=(fHistElectronCutLow[icut]->GetXaxis()->GetNbins()+1) &&
 494        etaBin!=0  && etaBin!=(fHistElectronCutLow[icut]->GetYaxis()->GetNbins()+1))
 495       lowElectronCut = fHistElectronCutLow[icut]->GetBinContent(centBin, etaBin, pbin);
 496   }
 497   if((fPartType[icut]==AliPID::kElectron) && fHistElectronCutUp[icut]) {
 498     Int_t pbin = fHistElectronCutUp[icut]->GetZaxis()->FindBin(mom);
 499     Int_t centBin = fHistElectronCutUp[icut]->GetXaxis()->FindBin(centrality);
 500     Int_t etaBin = fHistElectronCutUp[icut]->GetYaxis()->FindBin(eta);
 501     // apply this cut only in the p range existing in the histogram range
 502     if(pbin!=0    && pbin!=(fHistElectronCutUp[icut]->GetZaxis()->GetNbins()+1)             &&
 503        centBin!=0 && centBin!=(fHistElectronCutUp[icut]->GetXaxis()->GetNbins()+1) &&
 504        etaBin!=0  && etaBin!=(fHistElectronCutUp[icut]->GetYaxis()->GetNbins()+1))
 505       upElectronCut = fHistElectronCutUp[icut]->GetBinContent(centBin, etaBin, pbin);
 506   }
 507   Bool_t selected=((numberOfSigmas>=lowElectronCut)&&(numberOfSigmas<=upElectronCut))^fExclude[icut];
 508   return selected;
 509 }
 510
 511 //______________________________________________
 512 Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
 513 {
 514   //
 515   // TRD part of the pid check
 516   // the TRD checks on the probabilities.
 517   //
 518
 519   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 520   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 521   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 522
 523   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
 524
 525   Double_t p[AliPID::kSPECIES]={0.};
 526   fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
 527   Float_t particleProb=p[fPartType[icut]];
 528
 529   Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
 530   return selected;
 531 }
 532
 533 //______________________________________________
 534 Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod)
 535 {
 536   //
 537   // TRD part of the pid check using electron efficiency requirement
 538   // in this case the upper limit as well as the particle specie is ignored
 539   //   and the lower limit regarded as the requested electron efficiency
 540   //
 541
 542   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 543   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 544   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 545
 546   Double_t centrality = -1.;
 547   if(part->IsA() == AliESDtrack::Class())
 548     centrality=(const_cast<AliESDEvent*>( (static_cast<const AliESDtrack*>(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 549   if(part->IsA() == AliAODTrack::Class())
 550     centrality=(const_cast<AliAODEvent*>( (static_cast<const AliAODTrack*>(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 551
 552   Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut];
 553   return selected;
 554 }
 555
 556 //______________________________________________
 557 Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
 558 {
 559   //
 560   // TOF part of the PID check
 561   // Don't accept the track if there was no pid bit set
 562   //
 563   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,part);
 564   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 565   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 566
 567   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
 568
 569   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 570   return selected;
 571 }
 572
 573 //______________________________________________
 574 Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
 575 {
 576   //
 577   // emcal pid selecttion
 578   //
 579
 580   //TODO: correct way to check for emcal pid?
 581   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);
 582
 583   Bool_t hasPID=numberOfSigmas>-998.;
 584
 585   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
 586   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;
 587
 588
 589   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 590   return selected;
 591 }
 592
 593 //______________________________________________
 594 void AliDielectronPID::SetDefaults(Int_t def){
 595   //
 596   // initialise default pid strategies
 597   //
 598
 599   if (def==0){
 600     // 2sigma bands TPC:
 601     // - include e
 602     // - exclude mu,K,pi,p
 603     // -complete p range
 604     AddCut(kTPC,AliPID::kElectron,2);
 605     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
 606     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
 607     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
 608     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
 609   } else if (def==1) {
 610     // 2sigma bands TPC:
 611     // - include e 0<p<inf
 612     // - exclude mu,K,pi,p 0<p<2
 613     AddCut(kTPC,AliPID::kElectron,2.);
 614     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
 615     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
 616     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
 617     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
 618
 619   } else if (def==2) {
 620     // include 2sigma e TPC
 621     // 3sigma bands TOF
 622     // - exclude K,P
 623     AddCut(kTPC,AliPID::kElectron,-3.,3.);
 624     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 625     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 626
 627   } else if (def==3 || def==4) { // def3 and def4 are the same??
 628     // include 2sigma e TPC
 629     // 3sigma bands TOF
 630     // - exclude K 0<p<1
 631     // - exclude P 0<p<2
 632     AddCut(kTPC,AliPID::kElectron,2);
 633     AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
 634     AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
 635     AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
 636
 637   } else if (def==5) {
 638     AddCut(kTPC,AliPID::kElectron,-0.5,3);
 639     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 640   } else if (def==6) {
 641     // lower cut TPC: parametrisation by HFE
 642     // upper cut TPC: 3 sigma
 643     // TOF ele band 3sigma 0<p<1.5GeV
 644     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 645     lowerCut->SetParameters(-2.7,-0.4357);
 646     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 647     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 648   } else if (def==7) {
 649     // wide TPC cut
 650     // TOF ele band 3sigma 0<p<1.5GeV
 651     AddCut(kTPC,AliPID::kElectron,10.);
 652     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 653   } else if (def==8) {
 654     // TOF 5 sigma inclusion if TOFpid available
 655     // this should reduce K,p,Pi to a large extent
 656     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 657   } else if (def==9) {
 658     // lower cut TPC: parametrisation by HFE
 659     // upper cut TPC: 3 sigma
 660     // TOF 5 sigma inclusion if TOFpid available
 661     // this should reduce K,p,Pi to a large extent
 662     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 663     lowerCut->SetParameters(-2.65,-0.6757);
 664     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 665     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 666
 667   } else if (def==10) {
 668     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 669     AddCut(kTPC,AliPID::kElectron,3.);
 670     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 671     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 672
 673   } else if (def==11) {
 674     // lower cut TPC: parametrisation by HFE
 675     // only use from period d on !!
 676     // upper cut TPC: 3 sigma
 677     // TOF ele band 3sigma 0<p<2.0GeV
 678     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 679     lowerCut->SetParameters(-3.718,-0.4,0.27);
 680     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 681     AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
 682   } else if (def==12) {
 683     // lower cut TPC: parametrisation by HFE
 684     // only use from period d on !!
 685     // upper cut TPC: 3 sigma
 686     // TOF 5 sigma inclusion if TOFpid available
 687     // this should reduce K,p,Pi to a large extent
 688     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 689     lowerCut->SetParameters(-3.718,-0.4,0.27);
 690     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 691     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 692   }
 693   else if (def==13) {
 694     // TPC electron inclusion
 695     // TOF electron inclusion if available
 696     AddCut(kTOF,AliPID::kElectron,-4.,4.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 697     AddCut(kTPC,AliPID::kElectron,-3.5,3.5);
 698   }
 699   else if (def==14) {
 700     // TRD 1D 90% elec eff, 4-6 tracklets
 701     // TPC electron inclusion
 702     // TOF electron inclusion if available
 703     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 704            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 705     AddCut(kTOF,AliPID::kElectron,-4.,4.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 706     AddCut(kTPC,AliPID::kElectron,-3.5,3.5);
 707   }
 708   else if (def==15) {
 709     // TRD 1D 90% elec eff, 4-6 tracklets, chi2 < 2
 710     // TPC electron inclusion
 711     // TOF electron inclusion if available
 712     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 713            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 714     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,0.,2.,kFALSE,
 715            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDchi2);
 716     AddCut(kTOF,AliPID::kElectron,-4.,4.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 717     AddCut(kTPC,AliPID::kElectron,-3.5,3.5);
 718   }
 719
 720 }
 721
 722 //______________________________________________
 723 void AliDielectronPID::SetCorrVal(Double_t run)
 724 {
 725   //
 726   // set correction value for run
 727   //
 728   fgCorr=0.;
 729   fgCorrdEdx=1.;
 730
 731   if (fgFitCorr){
 732     fgCorr=fgFitCorr->Eval(run);
 733     if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
 734     if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
 735   }
 736
 737   if (fgdEdxRunCorr){
 738     fgCorrdEdx=fgdEdxRunCorr->Eval(run);
 739     if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
 740     if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
 741   }
 742 }
 743
 744 //______________________________________________
 745 Double_t AliDielectronPID::GetEtaCorr(const AliVTrack *track)
 746 {
 747   //
 748   // return eta correction
 749   //
 750   if (!fgFunEtaCorr) return 1;
 751   return fgFunEtaCorr->Eval(track->Eta());
 752 }
 753
 754 //______________________________________________
 755 void AliDielectronPID::SetCentroidCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
 756 {
 757   fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
 758   fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
 759   fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
 760   fgFunCntrdCorr=fun;
 761 }
 762 //______________________________________________
 763 void AliDielectronPID::SetWidthCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
 764 {
 765   fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
 766   fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
 767   fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
 768   fgFunWdthCorr=fun;
 769 }
 770
 771 //______________________________________________
 772 Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TF1 *fun)
 773 {
 774   //
 775   // return correction value
 776   //
 777
 778   //Fill only event and vparticle values (otherwise we end up in a circle)
 779   Double_t values[AliDielectronVarManager::kNMaxValues];
 780   AliDielectronVarManager::FillVarVParticle(track,values);
 781
 782   Int_t dim=fun->GetNdim();
 783   Double_t var[3] = {0.,0.,0.};
 784   if(dim>0) var[0] = values[fun->GetHistogram()->GetXaxis()->GetUniqueID()];
 785   if(dim>1) var[1] = values[fun->GetHistogram()->GetYaxis()->GetUniqueID()];
 786   if(dim>2) var[2] = values[fun->GetHistogram()->GetZaxis()->GetUniqueID()];
 787   Double_t corr = fun->Eval(var[0],var[1],var[2]);
 788   // printf(" %d-dim CORR value: %f (track %p) \n",dim,corr,track);
 789   return corr;
 790 }