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