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