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 <AliLog.h>
  34 #include <AliExternalTrackParam.h>
  35 #include <AliESDtrack.h>
  36 #include <AliESDpid.h>
  37 #include <AliAODpidUtil.h>
  38 #include <AliAODPid.h>
  39
  40 #include "AliDielectronVarManager.h"
  41
  42 #include "AliDielectronPID.h"
  43
  44 ClassImp(AliDielectronPID)
  45
  46 TGraph *AliDielectronPID::fgFitCorr=0x0;
  47 Double_t AliDielectronPID::fgCorr=0.0;
  48
  49 AliDielectronPID::AliDielectronPID() :
  50   AliAnalysisCuts(),
  51   fNcuts(0),
  52   fESDpid(0x0),
  53   fAODpidUtil(0x0)
  54 {
  55   //
  56   // Default Constructor
  57   //
  58   for (Int_t icut=0; icut<kNmaxPID; ++icut){
  59     fDetType[icut]=kTPC;
  60     fPartType[icut]=AliPID::kPion;
  61     fNsigmaLow[icut]=0;
  62     fNsigmaUp[icut]=0;
  63     fPmin[icut]=0;
  64     fPmax[icut]=0;
  65     fExclude[icut]=kFALSE;
  66     fFunUpperCut[icut]=0x0;
  67     fFunLowerCut[icut]=0x0;
  68     fRequirePIDbit[icut]=0;
  69   }
  70 }
  71
  72 //______________________________________________
  73 AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
  74   AliAnalysisCuts(name, title),
  75   fNcuts(0),
  76   fESDpid(0x0),
  77   fAODpidUtil(0x0)
  78 {
  79   //
  80   // Named Constructor
  81   //
  82   for (Int_t icut=0; icut<kNmaxPID; ++icut){
  83     fDetType[icut]=kTPC;
  84     fPartType[icut]=AliPID::kPion;
  85     fNsigmaLow[icut]=0;
  86     fNsigmaUp[icut]=0;
  87     fPmin[icut]=0;
  88     fPmax[icut]=0;
  89     fExclude[icut]=kFALSE;
  90     fFunUpperCut[icut]=0x0;
  91     fFunLowerCut[icut]=0x0;
  92     fRequirePIDbit[icut]=0;
  93   }
  94 }
  95
  96 //______________________________________________
  97 AliDielectronPID::~AliDielectronPID()
  98 {
  99   //
 100   // Default Destructor
 101   //
 102 }
 103
 104 //______________________________________________
 105 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp/*=-99999.*/,
 106                               Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
 107                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 108 {
 109   //
 110   // Add a pid nsigma cut
 111   // use response of detector 'det' in the momentum range ('pMin') to ['pMax']
 112   // use a sigma band betwee 'nSigmaLow' and 'nSigmaUp'
 113   // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
 114   // specify whether to 'exclude' the given band
 115   //
 116
 117   if (fNcuts>=kNmaxPID){
 118     AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
 119     return;
 120   }
 121   if (TMath::Abs(nSigmaUp+99999.)<1e-20){
 122     nSigmaUp=TMath::Abs(nSigmaLow);
 123     nSigmaLow=-1.*nSigmaUp;
 124   }
 125   fDetType[fNcuts]=det;
 126   fPartType[fNcuts]=type;
 127   fNsigmaLow[fNcuts]=nSigmaLow;
 128   fNsigmaUp[fNcuts]=nSigmaUp;
 129   fPmin[fNcuts]=pMin;
 130   fPmax[fNcuts]=pMax;
 131   fExclude[fNcuts]=exclude;
 132   fRequirePIDbit[fNcuts]=pidBitType;
 133   ++fNcuts;
 134
 135 }
 136
 137 //______________________________________________
 138 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
 139                               Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
 140                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 141 {
 142   //
 143   // cut using a TF1 as upper cut
 144   //
 145   if (funUp==0x0){
 146     AliError("A valid function is required for the upper cut. Not adding the cut!");
 147     return;
 148   }
 149   fFunUpperCut[fNcuts]=funUp;
 150   AddCut(det,type,nSigmaLow,0.,pMin,pMax,exclude,pidBitType);
 151 }
 152
 153 //______________________________________________
 154 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
 155                               Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
 156                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 157 {
 158   //
 159   // cut using a TF1 as lower cut
 160   //
 161   if (funLow==0x0){
 162     AliError("A valid function is required for the lower cut. Not adding the cut!");
 163     return;
 164   }
 165   fFunLowerCut[fNcuts]=funLow;
 166   AddCut(det,type,0.,nSigmaUp,pMin,pMax,exclude,pidBitType);
 167 }
 168
 169 //______________________________________________
 170 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
 171                               Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/,
 172                               UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
 173 {
 174   //
 175   // cut using a TF1 as lower and upper cut
 176   //
 177   if ( (funUp==0x0) || (funLow==0x0) ){
 178     AliError("A valid function is required for upper and lower cut. Not adding the cut!");
 179     return;
 180   }
 181   fFunUpperCut[fNcuts]=funUp;
 182   fFunLowerCut[fNcuts]=funLow;
 183   AddCut(det,type,0.,0.,pMin,pMax,exclude,pidBitType);
 184 }
 185
 186 //______________________________________________
 187 Bool_t AliDielectronPID::IsSelected(TObject* track)
 188 {
 189   //
 190   // perform PID cuts
 191   //
 192
 193   //loop over all cuts
 194   AliVTrack *part=static_cast<AliVTrack*>(track);
 195   //TODO: Which momentum to use?
 196   //      Different momenta for different detectors?
 197   Double_t mom=part->P();
 198
 199   Bool_t selected=kFALSE;
 200   fESDpid=AliDielectronVarManager::GetESDpid();
 201   fAODpidUtil=AliDielectronVarManager::GetAODpidUtil();
 202
 203   for (UChar_t icut=0; icut<fNcuts; ++icut){
 204     Double_t pMin=fPmin[icut];
 205     Double_t pMax=fPmax[icut];
 206
 207     // test momentum range. In case pMin==pMax use all momenta
 208     if ( (TMath::Abs(pMin-pMax)>1e-20) && (mom<=pMin || mom>pMax) ) continue;
 209
 210
 211     switch (fDetType[icut]){
 212     case kITS:
 213       selected = IsSelectedITS(part,icut);
 214       break;
 215     case kTPC:
 216       selected = IsSelectedTPC(part,icut);
 217       break;
 218     case kTRD:
 219       selected = IsSelectedTRD(part,icut);
 220       break;
 221     case kTOF:
 222       selected = IsSelectedTOF(part,icut);
 223       break;
 224     }
 225     if (!selected) return kFALSE;
 226   }
 227
 228   return selected;
 229 }
 230
 231 //______________________________________________
 232 Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
 233 {
 234   //
 235   // ITS part of the PID check
 236   // Don't accept the track if there was no pid bit set
 237   //
 238   Float_t numberOfSigmas=-1000.;
 239
 240   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kFALSE;
 241   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kTRUE;
 242
 243   Double_t mom=part->P();
 244
 245   if (part->IsA()==AliESDtrack::Class()){
 246     // ESD case in case the PID bit is not set, don't use this track!
 247     AliESDtrack *track=static_cast<AliESDtrack*>(part);
 248     numberOfSigmas=fESDpid->NumberOfSigmasITS(track, fPartType[icut]);
 249   }else if(part->IsA()==AliAODTrack::Class()){
 250     // AOD case
 251     // FIXME: Is there a place to check whether the PID is was set in ESD???
 252     AliAODTrack *track=static_cast<AliAODTrack*>(part);
 253     numberOfSigmas=fAODpidUtil->NumberOfSigmasITS(track, fPartType[icut]);
 254   }
 255
 256   // test if we are supposed to use a function for the cut
 257   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
 258   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
 259
 260   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 261   return selected;
 262 }
 263
 264 //______________________________________________
 265 Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut)
 266 {
 267   //
 268   // TPC part of the PID check
 269   // Don't accept the track if there was no pid bit set
 270   //
 271   Float_t numberOfSigmas=-1000.;
 272
 273   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kFALSE;
 274   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kTRUE;
 275
 276   Double_t mom=part->P();
 277
 278   if (part->IsA()==AliESDtrack::Class()){
 279     // ESD case in case the PID bit is not set, don't use this track!
 280     AliESDtrack *track=static_cast<AliESDtrack*>(part);
 281     const AliExternalTrackParam *in = track->GetInnerParam();
 282     if (in) mom = in->GetP();
 283     numberOfSigmas=fESDpid->NumberOfSigmasTPC(track, fPartType[icut]);
 284   }else if(part->IsA()==AliAODTrack::Class()){
 285     // AOD case
 286     // FIXME: Is there a place to check whether the PID is was set in ESD???
 287     AliAODTrack *track=static_cast<AliAODTrack*>(part);
 288     const AliAODPid *pidObj = track->GetDetPid();
 289     if (pidObj) mom = pidObj->GetTPCmomentum();
 290     numberOfSigmas=fAODpidUtil->NumberOfSigmasTPC(track, fPartType[icut]);
 291   }
 292   if (fPartType[icut]==AliPID::kElectron){
 293     numberOfSigmas-=fgCorr;
 294   }
 295
 296   // test if we are supposed to use a function for the cut
 297   if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
 298   if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
 299
 300   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 301   return selected;
 302 }
 303
 304 //______________________________________________
 305 Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /*part*/, Int_t /*icut*/)
 306 {
 307   //
 308   // TRD part of the pid check
 309   //
 310   return kFALSE;
 311 }
 312
 313 //______________________________________________
 314 Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
 315 {
 316   //
 317   // TOF part of the PID check
 318   // Don't accept the track if there was no pid bit set
 319   //
 320   Float_t numberOfSigmas=-1000.;
 321
 322   if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kTOFpid)) return kFALSE;
 323   if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTOFpid)) return kTRUE;
 324
 325   if (part->IsA()==AliESDtrack::Class()){
 326     // ESD case in case the PID bit is not set, don't use this track!
 327     AliESDtrack *track=static_cast<AliESDtrack*>(part);
 328     numberOfSigmas=fESDpid->NumberOfSigmasTOF(track, fPartType[icut], fESDpid->GetTOFResponse().GetTimeZero());
 329   }else if(part->IsA()==AliAODTrack::Class()){
 330     // AOD case
 331     // FIXME: Is there a place to check whether the PID is was set in ESD???
 332     AliAODTrack *track=static_cast<AliAODTrack*>(part);
 333     numberOfSigmas=fAODpidUtil->NumberOfSigmasTOF(track, fPartType[icut]);
 334   }
 335
 336   Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
 337   return selected;
 338 }
 339
 340 //______________________________________________
 341 void AliDielectronPID::SetDefaults(Int_t def){
 342   //
 343   // initialise default pid strategies
 344   //
 345
 346   if (def==0){
 347     // 2sigma bands TPC:
 348     // - include e
 349     // - exclude mu,K,pi,p
 350     // -complete p range
 351     AddCut(kTPC,AliPID::kElectron,2);
 352     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
 353     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
 354     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
 355     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
 356   } else if (def==1) {
 357     // 2sigma bands TPC:
 358     // - include e 0<p<inf
 359     // - exclude mu,K,pi,p 0<p<2
 360
 361     AddCut(kTPC,AliPID::kElectron,2.);
 362     AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
 363     AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
 364     AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
 365     AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
 366
 367   } else if (def==2) {
 368     // include 2sigma e TPC
 369     // 3sigma bands TOF
 370     // - exclude K,P
 371     AddCut(kTPC,AliPID::kElectron,-3.,3.);
 372     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 373     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 374
 375   } else if (def==3) {
 376     // include 2sigma e TPC
 377     // 3sigma bands TOF
 378     // - exclude K 0<p<1
 379     // - exclude P 0<p<2
 380     AddCut(kTPC,AliPID::kElectron,2);
 381     AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
 382     AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
 383     AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
 384
 385   } else if (def==4) {
 386     // include 2sigma e TPC
 387     // 3sigma bands TOF
 388     // - exclude K 0<p<1
 389     // - exclude P 0<p<2
 390     AddCut(kTPC,AliPID::kElectron,2.);
 391     AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
 392     AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
 393     AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
 394   } else if (def==5) {
 395     AddCut(kTPC,AliPID::kElectron,-0.5,3);
 396     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 397   } else if (def==6) {
 398     // lower cut TPC: parametrisation by HFE
 399     // upper cut TPC: 3 sigma
 400     // TOF ele band 3sigma 0<p<1.5GeV
 401     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 402     lowerCut->SetParameters(-2.7,-0.4357);
 403     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 404     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 405   } else if (def==7) {
 406     // wide TPC cut
 407     // TOF ele band 3sigma 0<p<1.5GeV
 408     AddCut(kTPC,AliPID::kElectron,10.);
 409     AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
 410   } else if (def==8) {
 411     // TOF 5 sigma inclusion if TOFpid available
 412     // this should reduce K,p,Pi to a large extent
 413     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 414   } else if (def==9) {
 415     // lower cut TPC: parametrisation by HFE
 416     // upper cut TPC: 3 sigma
 417     // TOF 5 sigma inclusion if TOFpid available
 418     // this should reduce K,p,Pi to a large extent
 419     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
 420     lowerCut->SetParameters(-2.65,-0.6757);
 421     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 422     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 423
 424   } else if (def==10) {
 425     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 426     AddCut(kTPC,AliPID::kElectron,3.);
 427     AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
 428     AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
 429
 430   } else if (def==11) {
 431     // lower cut TPC: parametrisation by HFE
 432     // only use from period d on !!
 433     // upper cut TPC: 3 sigma
 434     // TOF ele band 3sigma 0<p<2.0GeV
 435     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 436     lowerCut->SetParameters(-3.718,-0.4,0.27);
 437     AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
 438     AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
 439   } else if (def==12) {
 440     // lower cut TPC: parametrisation by HFE
 441     // only use from period d on !!
 442     // upper cut TPC: 3 sigma
 443     // TOF 5 sigma inclusion if TOFpid available
 444     // this should reduce K,p,Pi to a large extent
 445     TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
 446     lowerCut->SetParameters(-3.718,-0.4,0.27);
 447     AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
 448     AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
 449   }
 450
 451 }
 452
 453 //______________________________________________
 454 void AliDielectronPID::SetCorrVal(Double_t run)
 455 {
 456   //
 457   // set correction value for run
 458   //
 459   fgCorr=0;
 460   if (!fgFitCorr) return;
 461   fgCorr=fgFitCorr->Eval(run);
 462   if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
 463   if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
 464 //   printf("Corr: %f\n",fgCorr);
 465 }
 466