1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////
26 ///////////////////////////////////////////////////////////////////////////
33 #include <AliVTrack.h>
34 #include <AliVCluster.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>
43 #include "AliDielectronVarManager.h"
44 #include "AliDielectronVarCuts.h"
46 #include "AliDielectronPID.h"
48 ClassImp(AliDielectronPID)
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;
58 AliDielectronPID::AliDielectronPID() :
60 fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
65 // Default Constructor
67 for (Int_t icut=0; icut<kNmaxPID; ++icut){
69 fPartType[icut]=AliPID::kPion;
74 fExclude[icut]=kFALSE;
75 fFunUpperCut[icut]=0x0;
76 fFunLowerCut[icut]=0x0;
77 fRequirePIDbit[icut]=0;
83 fMapElectronCutLow[icut]=0x0;
87 //______________________________________________
88 AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
89 AliAnalysisCuts(name, title),
90 fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
97 for (Int_t icut=0; icut<kNmaxPID; ++icut){
99 fPartType[icut]=AliPID::kPion;
104 fExclude[icut]=kFALSE;
105 fFunUpperCut[icut]=0x0;
106 fFunLowerCut[icut]=0x0;
107 fRequirePIDbit[icut]=0;
109 fSigmaFunLow[icut]=0;
113 fMapElectronCutLow[icut]=0x0;
117 //______________________________________________
118 AliDielectronPID::~AliDielectronPID()
121 // Default Destructor
123 if (fUsedVars) delete fUsedVars;
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*/)
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
139 if (fNcuts>=kNmaxPID){
140 AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
143 if (TMath::Abs(nSigmaUp+99999.)<1e-20){
144 nSigmaUp=TMath::Abs(nSigmaLow);
145 nSigmaLow=-1.*nSigmaUp;
147 fDetType[fNcuts]=det;
148 fPartType[fNcuts]=type;
149 fNsigmaLow[fNcuts]=nSigmaLow;
150 fNsigmaUp[fNcuts]=nSigmaUp;
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);
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])));
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*/)
172 // cut using a TF1 as upper cut
175 AliError("A valid function is required for the upper cut. Not adding the cut!");
178 fFunUpperCut[fNcuts]=funUp;
179 AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var);
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*/)
188 // cut using a TF1 as lower cut
191 AliError("A valid function is required for the lower cut. Not adding the cut!");
194 fFunLowerCut[fNcuts]=funLow;
195 AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
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*/)
204 // cut using a TF1 as lower and upper cut
206 if ( (funUp==0x0) || (funLow==0x0) ){
207 AliError("A valid function is required for upper and lower cut. Not adding the cut!");
210 fFunUpperCut[fNcuts]=funUp;
211 fFunLowerCut[fNcuts]=funLow;
212 AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var);
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)
221 // cut using a TF1 as lower cut
224 AliError("A valid function is required for the sigma cut. Not adding the cut!");
227 fFunSigma[fNcuts]=funSigma;
228 fSigmaFunLow[fNcuts]=min;
229 fSigmaFunUp[fNcuts]=max;
231 AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1);
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*/)
240 // cut using a THnBase as a lower cut
243 AliError("A valid histogram is required for the lower cut. Not adding the cut!");
246 fMapElectronCutLow[fNcuts]=histLow;
247 AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var);
249 //______________________________________________
250 void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp,
251 AliDielectronVarCuts *var, Bool_t exclude/*=kFALSE*/,
252 UInt_t pidBitType/*=AliDielectronPID::kRequire*/)
255 // Add a pid nsigma cut
256 // use response of detector 'det' in the ranges for variables defined in var
257 // use a sigma band between 'nSigmaLow' and 'nSigmaUp'
258 // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
259 // specify whether to 'exclude' the given band
262 if (fNcuts>=kNmaxPID){
263 AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
266 if (TMath::Abs(nSigmaUp+99999.)<1e-20){
267 nSigmaUp=TMath::Abs(nSigmaLow);
268 nSigmaLow=-1.*nSigmaUp;
270 fDetType[fNcuts]=det;
271 fPartType[fNcuts]=type;
272 fNsigmaLow[fNcuts]=nSigmaLow;
273 fNsigmaUp[fNcuts]=nSigmaUp;
274 fExclude[fNcuts]=exclude;
275 fRequirePIDbit[fNcuts]=pidBitType;
276 fVarCuts[fNcuts]=var;
278 AliDebug(1,Form("Add PID cut %d: sigma [% .1f,% .1f] \n",
279 fNcuts,nSigmaLow,nSigmaUp));
285 //______________________________________________
286 Bool_t AliDielectronPID::IsSelected(TObject* track)
293 AliVTrack *part=static_cast<AliVTrack*>(track);
294 AliESDtrack *esdTrack=0x0;
295 AliAODTrack *aodTrack=0x0;
296 Double_t origdEdx=-1;
298 // apply ETa correction, remove once this is in the tender
299 if( (part->IsA() == AliESDtrack::Class()) ){
300 esdTrack=static_cast<AliESDtrack*>(part);
301 origdEdx=esdTrack->GetTPCsignal();
302 esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
303 } else if ( (part->IsA() == AliAODTrack::Class()) ){
304 aodTrack=static_cast<AliAODTrack*>(track);
305 AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
307 origdEdx=pid->GetTPCsignal();
308 pid->SetTPCsignal(origdEdx/GetEtaCorr(aodTrack)/fgCorrdEdx);
312 // check for corrections and add their variables to the fill map
314 fUsedVars->SetBitNumber(fgFunCntrdCorr->GetXaxis()->GetUniqueID(), kTRUE);
315 fUsedVars->SetBitNumber(fgFunCntrdCorr->GetYaxis()->GetUniqueID(), kTRUE);
316 fUsedVars->SetBitNumber(fgFunCntrdCorr->GetZaxis()->GetUniqueID(), kTRUE);
319 fUsedVars->SetBitNumber(fgFunWdthCorr->GetXaxis()->GetUniqueID(), kTRUE);
320 fUsedVars->SetBitNumber(fgFunWdthCorr->GetYaxis()->GetUniqueID(), kTRUE);
321 fUsedVars->SetBitNumber(fgFunWdthCorr->GetZaxis()->GetUniqueID(), kTRUE);
323 for(UChar_t icut=0; icut<fNcuts; ++icut) {
324 if(!fMapElectronCutLow[icut]) continue;
325 for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
326 TString var(fMapElectronCutLow[icut]->GetAxis(idim)->GetName());
327 fUsedVars->SetBitNumber(AliDielectronVarManager::GetValueType(var.Data()), kTRUE);
332 Double_t values[AliDielectronVarManager::kNMaxValues];
333 AliDielectronVarManager::SetFillMap(fUsedVars);
334 AliDielectronVarManager::Fill(track,values);
336 Bool_t selected=kFALSE;
337 fPIDResponse=AliDielectronVarManager::GetPIDResponse();
338 for (UChar_t icut=0; icut<fNcuts; ++icut){
339 Double_t min=fmin[icut];
340 Double_t max=fmax[icut];
341 Double_t val=values[fActiveCuts[icut]];
343 // test var range. In case min==max do not cut
344 if ( fVarCuts[icut] ) {
345 if ( !fVarCuts[icut]->IsSelected(part) ) continue;
347 else if ( ( (TMath::Abs(min-max)>1e-20) && (val<min || val>=max) ) ) {
351 // check if fFunSigma is set, then check if 'part' is in sigma range of the function
353 val= fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
354 if (fPartType[icut]==AliPID::kElectron){
357 min= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunLow[icut];
358 max= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunUp[icut];
359 if(val<min || val>=max) continue;
362 switch (fDetType[icut]){
364 selected = IsSelectedITS(part,icut);
367 selected = IsSelectedTPC(part,icut,values);
370 selected = IsSelectedTRD(part,icut);
373 selected = IsSelectedTRDeleEff(part,icut);
376 selected = IsSelectedTRDeleEff(part,icut,AliTRDPIDResponse::kLQ2D);
379 selected = IsSelectedTOF(part,icut);
382 selected = IsSelectedEMCAL(part,icut);
386 if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
388 AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
389 if (pid) pid->SetTPCsignal(origdEdx);
396 if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
398 AliAODPid *pid=const_cast<AliAODPid*>(aodTrack->GetDetPid());
399 if (pid) pid->SetTPCsignal(origdEdx);
404 //______________________________________________
405 Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut)
408 // ITS part of the PID check
409 // Don't accept the track if there was no pid bit set
411 AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kITS,part);
412 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
413 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
415 Double_t mom=part->P();
417 Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]);
419 // test if we are supposed to use a function for the cut
420 if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom);
421 if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom);
423 Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
427 //______________________________________________
428 Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut, Double_t *values)
431 // TPC part of the PID check
432 // Don't accept the track if there was no pid bit set
434 AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC,part);
435 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
436 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
439 Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]);
442 if (fPartType[icut]==AliPID::kElectron){
444 numberOfSigmas-=fgCorr;
445 // via functions (1-3D)
446 numberOfSigmas-=GetCntrdCorr(part);
447 numberOfSigmas/=GetWdthCorr(part);
450 // matching of MC and data //
452 // test if we are supposed to use a function for the cut (matching via fits)
453 if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
454 if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(values[AliDielectronVarManager::kPIn]);
456 // test if we are using cut maps (in calibrated units)
457 Double_t lowElectronCut = fNsigmaLow[icut];
458 Double_t upElectronCut = fNsigmaUp[icut];
460 if((fPartType[icut]==AliPID::kElectron) && (fMapElectronCutLow[icut] )) {
461 Double_t *vals = new Double_t[fMapElectronCutLow[icut]->GetNdimensions()];//={-1};
462 // get array of values for the corresponding dimensions using axis names
463 for(Int_t idim=0; idim<fMapElectronCutLow[icut]->GetNdimensions(); idim++) {
464 vals[idim] = values[AliDielectronVarManager::GetValueType(fMapElectronCutLow[icut]->GetAxis(idim)->GetName())];
465 // printf(" \t %s %.3f ",fMapElectronCutLow[icut]->GetAxis(idim)->GetName(),vals[idim]);
467 // find bin for values (w/o creating it in case it is not filled)
468 Long_t bin = fMapElectronCutLow[icut]->GetBin(vals,kFALSE);
469 if(bin>0) lowElectronCut=fMapElectronCutLow[icut]->GetBinContent(bin);
470 else lowElectronCut=100;
471 // printf(" low cut %.3f (%ld) \n", lowElectronCut,bin);
476 Bool_t selected=((numberOfSigmas>=lowElectronCut)&&(numberOfSigmas<=upElectronCut))^fExclude[icut];
480 //______________________________________________
481 Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut)
484 // TRD part of the pid check
485 // the TRD checks on the probabilities.
488 AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
489 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
490 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
492 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE;
494 Double_t p[AliPID::kSPECIES]={0.};
495 fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p);
496 Float_t particleProb=p[fPartType[icut]];
498 Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut];
502 //______________________________________________
503 Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod)
506 // TRD part of the pid check using electron efficiency requirement
507 // in this case the upper limit as well as the particle specie is ignored
508 // and the lower limit regarded as the requested electron efficiency
511 AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTRD,part);
512 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
513 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
515 Double_t centrality = -1.;
516 if(part->IsA() == AliESDtrack::Class())
517 centrality=(const_cast<AliESDEvent*>( (static_cast<const AliESDtrack*>(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
518 if(part->IsA() == AliAODTrack::Class())
519 centrality=(const_cast<AliAODEvent*>( (static_cast<const AliAODTrack*>(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M");
521 Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut];
525 //______________________________________________
526 Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut)
529 // TOF part of the PID check
530 // Don't accept the track if there was no pid bit set
532 AliPIDResponse::EDetPidStatus pidStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,part);
533 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kFALSE;
534 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&(pidStatus!=AliPIDResponse::kDetPidOk)) return kTRUE;
536 Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]);
538 Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
542 //______________________________________________
543 Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut)
546 // emcal pid selecttion
549 //TODO: correct way to check for emcal pid?
550 Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]);
552 Bool_t hasPID=numberOfSigmas>-998.;
554 if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE;
555 if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE;
558 Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut];
562 //______________________________________________
563 void AliDielectronPID::SetDefaults(Int_t def){
565 // initialise default pid strategies
571 // - exclude mu,K,pi,p
573 AddCut(kTPC,AliPID::kElectron,2);
574 AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE);
575 AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE);
576 AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE);
577 AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE);
580 // - include e 0<p<inf
581 // - exclude mu,K,pi,p 0<p<2
582 AddCut(kTPC,AliPID::kElectron,2.);
583 AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
584 AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
585 AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
586 AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
589 // include 2sigma e TPC
592 AddCut(kTPC,AliPID::kElectron,-3.,3.);
593 AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
594 AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
596 } else if (def==3 || def==4) { // def3 and def4 are the same??
597 // include 2sigma e TPC
601 AddCut(kTPC,AliPID::kElectron,2);
602 AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
603 AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
604 AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
607 AddCut(kTPC,AliPID::kElectron,-0.5,3);
608 AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
610 // lower cut TPC: parametrisation by HFE
611 // upper cut TPC: 3 sigma
612 // TOF ele band 3sigma 0<p<1.5GeV
613 TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
614 lowerCut->SetParameters(-2.7,-0.4357);
615 AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
616 AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
619 // TOF ele band 3sigma 0<p<1.5GeV
620 AddCut(kTPC,AliPID::kElectron,10.);
621 AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
623 // TOF 5 sigma inclusion if TOFpid available
624 // this should reduce K,p,Pi to a large extent
625 AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
627 // lower cut TPC: parametrisation by HFE
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)", 0, 100);
632 lowerCut->SetParameters(-2.65,-0.6757);
633 AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
634 AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
636 } else if (def==10) {
637 AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
638 AddCut(kTPC,AliPID::kElectron,3.);
639 AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE);
640 AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
642 } else if (def==11) {
643 // lower cut TPC: parametrisation by HFE
644 // only use from period d on !!
645 // upper cut TPC: 3 sigma
646 // TOF ele band 3sigma 0<p<2.0GeV
647 TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
648 lowerCut->SetParameters(-3.718,-0.4,0.27);
649 AddCut(kTPC,AliPID::kElectron,lowerCut,3.);
650 AddCut(kTOF,AliPID::kElectron,-3,3,0,2.);
651 } else if (def==12) {
652 // lower cut TPC: parametrisation by HFE
653 // only use from period d on !!
654 // upper cut TPC: 3 sigma
655 // TOF 5 sigma inclusion if TOFpid available
656 // this should reduce K,p,Pi to a large extent
657 TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100);
658 lowerCut->SetParameters(-3.718,-0.4,0.27);
659 AddCut(kTPC,AliPID::kElectron,lowerCut,4.);
660 AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable);
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,10.);
669 // TRD 1D 90% elec eff, 4-6 tracklets
670 // TPC electron inclusion
671 // TOF electron inclusion if available
672 AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
673 AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
674 AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
675 AddCut(kTPC,AliPID::kElectron,10.);
678 // TRD 1D 90% elec eff, 4-6 tracklets, chi2 < 2
679 // TPC electron inclusion
680 // TOF electron inclusion if available
681 AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,3.5,6.,kFALSE,
682 AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDpidQuality);
683 AddCut(AliDielectronPID::kTRDeleEff,AliPID::kElectron,.9,1.,0.,2.,kFALSE,
684 AliDielectronPID::kIfAvailable,AliDielectronVarManager::kTRDchi2);
685 AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
686 AddCut(kTPC,AliPID::kElectron,10.);
689 // TPC electron inclusion
690 // TOF electron inclusion if available
691 AddCut(kTOF,AliPID::kElectron,-3.,3.,0,200,kFALSE,AliDielectronPID::kIfAvailable);
692 AddCut(kTPC,AliPID::kElectron,-3.5,+3.5);
697 //______________________________________________
698 void AliDielectronPID::SetCorrVal(Double_t run)
701 // set correction value for run
707 fgCorr=fgFitCorr->Eval(run);
708 if (run<fgFitCorr->GetX()[0]) fgCorr=fgFitCorr->GetY()[0];
709 if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1];
713 fgCorrdEdx=fgdEdxRunCorr->Eval(run);
714 if (run<fgdEdxRunCorr->GetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0];
715 if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1];
719 //______________________________________________
720 Double_t AliDielectronPID::GetEtaCorr(const AliVTrack *track)
723 // return eta correction
725 if (!fgFunEtaCorr) return 1;
726 return fgFunEtaCorr->Eval(track->Eta());
729 //______________________________________________
730 Double_t AliDielectronPID::GetPIDCorr(const AliVTrack *track, TH1 *hist)
733 // return correction value
735 //TODO: think about taking an values array as argument to reduce # var fills
737 //Fill only event and vparticle values (otherwise we end up in a circle)
738 Double_t values[AliDielectronVarManager::kNMaxValues];
739 AliDielectronVarManager::FillVarVParticle(track,values);
741 TF1 *fun = (TF1*)hist->GetListOfFunctions()->At(0);
742 Int_t dim=fun->GetNdim();
744 Double_t var[3] = {0.,0.,0.};
745 if(dim>0) var[0] = values[hist->GetXaxis()->GetUniqueID()];
746 if(dim>1) var[1] = values[hist->GetYaxis()->GetUniqueID()];
747 if(dim>2) var[2] = values[hist->GetZaxis()->GetUniqueID()];
748 Double_t corr = fun->Eval(var[0],var[1],var[2]);
749 // printf("%d-dim CORR value: %f (track %p) \n",dim,corr,track);