major dielectron update (included also the data and plotting macros for paper)
[u/mrichter/AliRoot.git] / PWG3 / 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