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