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   // fill used variables into map
 248   for(Int_t idim=0; idim<fMapElectronCutLow[fNcuts]->GetNdimensions(); idim++) {
 249     TString vari(fMapElectronCutLow[fNcuts]->GetAxis(idim)->GetName());
 250     fUsedVars->SetBitNumber(AliDielectronVarManager::GetValueType(vari.Data()), kTRUE);
 251   }
 252   AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
 253 }
 254 //______________________________________________
 255 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
 256                               AliDielectronVarCuts *var, Bool_t exclude/*=kFALSE*/,
 257                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 258 {
 259   //
 260   // Add a pid nsigma cut
 261   // use response of detector 'det' in the ranges for variables defined in var
 262   // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
 263   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 264   // specify whether to 'exclude' the given band
 265   //
 266   if(!var) return;
 267   if (fNcuts>=kNmaxPID){
 268     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 269     return;
 270   }
 271   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 272     nSigmaUp=TMath::Abs(nSigmaLow);
 273     nSigmaLow=-1.*nSigmaUp;
 274   }
 275   fDetType[fNcuts]=det;
 276   fPartType[fNcuts]=type;
 277   fNsigmaLow[fNcuts]=nSigmaLow;
 278   fNsigmaUp[fNcuts]=nSigmaUp;
 279   fExclude[fNcuts]=exclude;
 280   fRequirePIDbit[fNcuts]=pidBitType;
 281   fVarCuts[fNcuts]=var;
 282
 283   AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \n",
 284                   fNcuts,nSigmaLow,nSigmaUp));
 285
 286   ++fNcuts;
 287
 288 }
 289
 290 //______________________________________________
 291 Bool_t AliDielectronPID::IsSelected(TObject* track)
 292 {
 293   //
 294   // perform PID cuts
 295   //
 296
 297   //loop over all cuts
 298   AliVTrack *part=static_cast<AliVTrack*>(track);
 299   AliESDtrack *esdTrack=0x0;
 300   AliAODTrack *aodTrack=0x0;
 301   Double_t origdEdx=-1;
 302
 303   // apply ETa correction, remove once this is in the tender
 304   if( (part->IsA() == AliESDtrack::Class()) ){
 305     esdTrack=static_cast<AliESDtrack*>(part);
 306     origdEdx=esdTrack->GetTPCsignal();
 307     esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 308   } else if ( (part->IsA() == AliAODTrack::Class()) ){
 309     aodTrack=static_cast<AliAODTrack*>(track);
 310     AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 311     if (pid){
 312       origdEdx=pid->GetTPCsignal();
 313       pid->SetTPCsignal(origdEdx/GetEtaCorr(aodTrack)/fgCorrdEdx);
 314     }
 315   }
 316
 317   // check for corrections and add their variables to the fill map
 318   if(fgFunCntrdCorr)  {
 319     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetXaxis()->GetUniqueID(), kTRUE);
 320     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetYaxis()->GetUniqueID(), kTRUE);
 321     fUsedVars->SetBitNumber(fgFunCntrdCorr->GetZaxis()->GetUniqueID(), kTRUE);
 322   }
 323   if(fgFunWdthCorr)  {
 324     fUsedVars->SetBitNumber(fgFunWdthCorr->GetXaxis()->GetUniqueID(), kTRUE);
 325     fUsedVars->SetBitNumber(fgFunWdthCorr->GetYaxis()->GetUniqueID(), kTRUE);
 326     fUsedVars->SetBitNumber(fgFunWdthCorr->GetZaxis()->GetUniqueID(), kTRUE);
 327   }
 328
 329   //Fill values
 330   Double_t values[AliDielectronVarManager::kNMaxValues];
 331   AliDielectronVarManager::SetFillMap(fUsedVars);
 332   AliDielectronVarManager::Fill(track,values);
 333
 334   Bool_t selected=kFALSE;
 335   fPIDResponse=AliDielectronVarManager::GetPIDResponse();
 336   for (UChar_t icut=0; icut<fNcuts; ++icut){
 337     Double_t min=fmin[icut];
 338     Double_t max=fmax[icut];
 339     Double_t val=values[fActiveCuts[icut]];
 340
 341     // test var range. In case min==max do not cut
 342     if ( fVarCuts[icut] ) {
 343       if ( !fVarCuts[icut]->IsSelected(part) ) continue;
 344     }
 345     else if ( ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) ) {
 346         continue;
 347     }
 348
 349     // check if fFunSigma is set, then check if 'part' is in sigma range of the function
 350     if(fFunSigma[icut]){
 351         val= fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
 352         if (fPartType[icut]==AliPID::kElectron){
 353             val-=fgCorr;
 354         }
 355         min= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunLow[icut];
 356         max= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunUp[icut];
 357         if(val<min || val>=max) continue;
 358     }
 359
 360     switch (fDetType[icut]){
 361     case kITS:
 362       selected = IsSelectedITS(part,icut);
 363       break;
 364     case kTPC:
 365       selected = IsSelectedTPC(part,icut,values);
 366       break;
 367     case kTRD:
 368       selected = IsSelectedTRD(part,icut);
 369       break;
 370     case kTRDeleEff:
 371       selected = IsSelectedTRDeleEff(part,icut);
 372       break;
 373     case kTRDeleEff2D:
 374       selected = IsSelectedTRDeleEff(part,icut,AliTRDPIDResponse::kLQ2D);
 375       break;
 376     case kTOF:
 377       selected = IsSelectedTOF(part,icut);
 378       break;
 379     case kEMCAL:
 380       selected = IsSelectedEMCAL(part,icut);
 381       break;
 382     }
 383     if (!selected) {
 384       if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 385       else if (aodTrack){
 386         AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 387         if (pid) pid->SetTPCsignal(origdEdx);
 388       }
 389
 390       return kFALSE;
 391     }
 392   }
 393
 394   if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 395   else if (aodTrack){
 396     AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
 397     if (pid) pid->SetTPCsignal(origdEdx);
 398   }
 399   return selected;
 400 }
 401
 402 //______________________________________________
 403 Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
 404 {
 405   //
 406   // ITS part of the PID check
 407   // Don't accept the track if there was no pid bit set
 408   //
 409   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kITS,part);
 410   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 411   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 412
 413   Double_t mom=part->P();
 414
 415   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]);
 416
 417   // test if we are supposed to use a function for the cut
 418   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
 419   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
 420
 421   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 422   return selected;
 423 }
 424
 425 //______________________________________________
 426 Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut, Double_t *values)
 427 {
 428   //
 429   // TPC part of the PID check
 430   // Don't accept the track if there was no pid bit set
 431   //
 432   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,part);
 433   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 434   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 435
 436
 437   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
 438
 439   // eta corrections
 440   if (fPartType[icut]==AliPID::kElectron){
 441     // old way 1D
 442     numberOfSigmas-=fgCorr;
 443     // via functions (1-3D)
 444     numberOfSigmas-=GetCntrdCorr(part);
 445     numberOfSigmas/=GetWdthCorr(part);
 446   }
 447
 448   // matching of MC and data //
 449
 450   // test if we are supposed to use a function for the cut (matching via fits)
 451   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
 452   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
 453
 454   // test if we are using cut maps (in calibrated units)
 455   Double_t lowElectronCut = fNsigmaLow[icut];
 456   Double_t upElectronCut  = fNsigmaUp[icut];
 457
 458   if((fPartType[icut]==AliPID::kElectron) && (fMapElectronCutLow[icut] )) {
 459     Double_t *vals = new Double_t[fMapElectronCutLow[icut]->GetNdimensions()];//={-1};
 460     // get array of values for the corresponding dimensions using axis names
 461     for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
 462       vals[idim] = values[AliDielectronVarManager::GetValueType(fMapElectronCutLow[icut]->GetAxis(idim)->GetName())];
 463       // printf(" \t %s %.3f ",fMapElectronCutLow[icut]->GetAxis(idim)->GetName(),vals[idim]);
 464     }
 465     // find bin for values (w/o creating it in case it is not filled)
 466     Long_t bin = fMapElectronCutLow[icut]->GetBin(vals,kFALSE);
 467     if(bin>0) lowElectronCut=fMapElectronCutLow[icut]->GetBinContent(bin);
 468     else lowElectronCut=100;
 469     //    printf("  low cut %.3f (%ld) \n", lowElectronCut,bin);
 470     delete [] vals;
 471   }
 472
 473
 474   Bool_t selected=((numberOfSigmas>=lowElectronCut)&&(numberOfSigmas<=upElectronCut))^fExclude[icut];
 475   return selected;
 476 }
 477
 478 //______________________________________________
 479 Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
 480 {
 481   //
 482   // TRD part of the pid check
 483   // the TRD checks on the probabilities.
 484   //
 485
 486   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 487   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 488   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 489
 490   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
 491
 492   Double_t p[AliPID::kSPECIES]={0.};
 493   fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
 494   Float_t particleProb=p[fPartType[icut]];
 495
 496   Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
 497   return selected;
 498 }
 499
 500 //______________________________________________
 501 Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod)
 502 {
 503   //
 504   // TRD part of the pid check using electron efficiency requirement
 505   // in this case the upper limit as well as the particle specie is ignored
 506   //   and the lower limit regarded as the requested electron efficiency
 507   //
 508
 509   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
 510   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 511   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 512
 513   Double_t centrality = -1.;
 514   if(part->IsA() == AliESDtrack::Class())
 515     centrality=(const_cast<AliESDEvent*>( (static_cast<const AliESDtrack*>(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 516   if(part->IsA() == AliAODTrack::Class())
 517     centrality=(const_cast<AliAODEvent*>( (static_cast<const AliAODTrack*>(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
 518
 519   Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut];
 520   return selected;
 521 }
 522
 523 //______________________________________________
 524 Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
 525 {
 526   //
 527   // TOF part of the PID check
 528   // Don't accept the track if there was no pid bit set
 529   //
 530   AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,part);
 531   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
 532   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
 533
 534   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
 535
 536   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 537   return selected;
 538 }
 539
 540 //______________________________________________
 541 Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
 542 {
 543   //
 544   // emcal pid selecttion
 545   //
 546
 547   //TODO: correct way to check for emcal pid?
 548   Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);
 549
 550   Bool_t hasPID=numberOfSigmas>-998.;
 551
 552   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
 553   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;
 554
 555
 556   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 557   return selected;
 558 }
 559
 560 //______________________________________________
 561 void AliDielectronPID::SetDefaults(Int_t def){
 562   //
 563   // initialise default pid strategies
 564   //
 565
 566   if (def==0){
 567     // 2sigma bands TPC:
 568     // - include e
 569     // - exclude mu,K,pi,p
 570     // -complete p range
 571     AddCut(kTPC,AliPID::kElectron,2);
 572     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
 573     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
 574     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
 575     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
 576   } else if (def==1) {
 577     // 2sigma bands TPC:
 578     // - include e 0<p<inf
 579     // - exclude mu,K,pi,p 0<p<2
 580     AddCut(kTPC,AliPID::kElectron,2.);
 581     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
 582     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
 583     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
 584     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
 585
 586   } else if (def==2) {
 587     // include 2sigma e TPC
 588     // 3sigma bands TOF
 589     // - exclude K,P
 590     AddCut(kTPC,AliPID::kElectron,-3.,3.);
 591     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 592     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 593
 594   } else if (def==3 || def==4) { // def3 and def4 are the same??
 595     // include 2sigma e TPC
 596     // 3sigma bands TOF
 597     // - exclude K 0<p<1
 598     // - exclude P 0<p<2
 599     AddCut(kTPC,AliPID::kElectron,2);
 600     AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
 601     AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
 602     AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
 603
 604   } else if (def==5) {
 605     AddCut(kTPC,AliPID::kElectron,-0.5,3);
 606     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 607   } else if (def==6) {
 608     // lower cut TPC: parametrisation by HFE
 609     // upper cut TPC: 3 sigma
 610     // TOF ele band 3sigma 0<p<1.5GeV
 611     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 612     lowerCut->SetParameters(-2.7,-0.4357);
 613     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 614     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 615   } else if (def==7) {
 616     // wide TPC cut
 617     // TOF ele band 3sigma 0<p<1.5GeV
 618     AddCut(kTPC,AliPID::kElectron,10.);
 619     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 620   } else if (def==8) {
 621     // TOF 5 sigma inclusion if TOFpid available
 622     // this should reduce K,p,Pi to a large extent
 623     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 624   } else if (def==9) {
 625     // lower cut TPC: parametrisation by HFE
 626     // upper cut TPC: 3 sigma
 627     // TOF 5 sigma inclusion if TOFpid available
 628     // this should reduce K,p,Pi to a large extent
 629     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 630     lowerCut->SetParameters(-2.65,-0.6757);
 631     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 632     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 633
 634   } else if (def==10) {
 635     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 636     AddCut(kTPC,AliPID::kElectron,3.);
 637     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 638     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 639
 640   } else if (def==11) {
 641     // lower cut TPC: parametrisation by HFE
 642     // only use from period d on !!
 643     // upper cut TPC: 3 sigma
 644     // TOF ele band 3sigma 0<p<2.0GeV
 645     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 646     lowerCut->SetParameters(-3.718,-0.4,0.27);
 647     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 648     AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
 649   } else if (def==12) {
 650     // lower cut TPC: parametrisation by HFE
 651     // only use from period d on !!
 652     // upper cut TPC: 3 sigma
 653     // TOF 5 sigma inclusion if TOFpid available
 654     // this should reduce K,p,Pi to a large extent
 655     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 656     lowerCut->SetParameters(-3.718,-0.4,0.27);
 657     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 658     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 659   }
 660   else if (def==13) {
 661     // TPC electron inclusion
 662     // TOF electron inclusion if available
 663     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 664     AddCut(kTPC,AliPID::kElectron,10.);
 665   }
 666   else if (def==14) {
 667     // TRD 1D 90% elec eff, 4-6 tracklets
 668     // TPC electron inclusion
 669     // TOF electron inclusion if available
 670     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 671            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 672     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 673     AddCut(kTPC,AliPID::kElectron,10.);
 674   }
 675   else if (def==15) {
 676     // TRD 1D 90% elec eff, 4-6 tracklets, chi2 < 2
 677     // TPC electron inclusion
 678     // TOF electron inclusion if available
 679     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
 680            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
 681     AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,0.,2.,kFALSE,
 682            AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDchi2);
 683     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 684     AddCut(kTPC,AliPID::kElectron,10.);
 685   }
 686   else if (def==16) {
 687     // TPC electron inclusion
 688     // TOF electron inclusion if available
 689     AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 690     AddCut(kTPC,AliPID::kElectron,-3.5,+3.5);
 691   }
 692
 693 }
 694
 695 //______________________________________________
 696 void AliDielectronPID::SetCorrVal(Double_t run)
 697 {
 698   //
 699   // set correction value for run
 700   //
 701   fgCorr=0.;
 702   fgCorrdEdx=1.;
 703
 704   if (fgFitCorr){
 705     fgCorr=fgFitCorr->Eval(run);
 706     if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
 707     if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
 708   }
 709
 710   if (fgdEdxRunCorr){
 711     fgCorrdEdx=fgdEdxRunCorr->Eval(run);
 712     if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
 713     if (run>fgdEdxRunCorr->GetX()[fgdEdxRunCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
 714   }
 715 }
 716
 717 //______________________________________________
 718 Double_t AliDielectronPID::GetEtaCorr(const AliVTrack *track)
 719 {
 720   //
 721   // return eta correction
 722   //
 723   if (!fgFunEtaCorr) return 1;
 724   return fgFunEtaCorr->Eval(track->Eta());
 725 }
 726
 727 //______________________________________________
 728 Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TH1 *hist)
 729 {
 730   //
 731   // return correction value
 732   //
 733   //TODO: think about taking an values array as argument to reduce # var fills
 734
 735   //Fill only event and vparticle values (otherwise we end up in a circle)
 736   Double_t values[AliDielectronVarManager::kNMaxValues];
 737   AliDielectronVarManager::FillVarVParticle(track,values);
 738
 739   TF1 *fun = (TF1*)hist->GetListOfFunctions()->At(0);
 740   Int_t dim=(fun?fun->GetNdim():hist->GetDimension());
 741
 742   Double_t var[3] = {0.,0.,0.};
 743   if(dim>0) var[0] = values[hist->GetXaxis()->GetUniqueID()];
 744   if(dim>1) var[1] = values[hist->GetYaxis()->GetUniqueID()];
 745   if(dim>2) var[2] = values[hist->GetZaxis()->GetUniqueID()];
 746   Double_t corr = 0.0;
 747   if(fun) corr = fun->Eval(var[0],var[1],var[2]);
 748   else    corr = hist->GetBinContent( hist->FindFixBin(var[0],var[1],var[2]) );
 749   //  for(Int_t i=0;i<dim;i++) printf("%d:%.3f  ",i,var[i]);
 750   //  printf("%d-dim CORR value: %f (track %p) \n",dim,corr,track);
 751   return corr;
 752 }