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