985884b65e3011b660c1f6f65956c206a7beedd2
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / AntiprotonToProton / AliProtonAnalysisBase.cxx
1 /**************************************************************************
2  * Author: Panos Christakoglou.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 /* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
15
16 //-----------------------------------------------------------------
17 //                 AliProtonAnalysisBase class
18 //   This is the class to deal with the proton analysis
19 //   Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21 #include <Riostream.h>
22 #include <TCanvas.h>
23 #include <TLatex.h>
24 #include <TF1.h>
25 #include <TList.h>
26 #include <TH1F.h>
27
28 #include <AliExternalTrackParam.h>
29 #include <AliESDEvent.h>
30 #include <AliPID.h>
31 #include <AliVertexerTracks.h>
32 #include <AliESDpid.h>
33 #include <AliTPCPIDResponse.h>
34 class AliLog;
35 class AliESDVertex;
36
37 #include "AliProtonAnalysisBase.h"
38
39 ClassImp(AliProtonAnalysisBase)
40
41 //____________________________________________________________________//
42 AliProtonAnalysisBase::AliProtonAnalysisBase() : 
43   TObject(),  fProtonAnalysisLevel("ESD"), fAnalysisMC(kFALSE),
44   fTriggerMode(kMB2), kUseOnlineTrigger(kFALSE), kUseOfflineTrigger(kFALSE), 
45   fPhysicsSelection(0),
46   fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
47   fAnalysisEtaMode(kFALSE),
48   fRunQAAnalysis(kFALSE),
49   fVxMax(100.), fVyMax(100.), fVzMax(100.), fMinNumOfContributors(0),
50   fNBinsX(0), fMinX(0.), fMaxX(0.),
51   fNBinsY(0), fMinY(0.), fMaxY(0.),
52   fMinTPCClusters(0), fMinITSClusters(0),
53   fMaxChi2PerTPCCluster(0.), fMaxChi2PerITSCluster(0.),
54   fMaxCov11(0.), fMaxCov22(0.), fMaxCov33(0.), fMaxCov44(0.), fMaxCov55(0.),
55   fMaxSigmaToVertex(0.), fMaxSigmaToVertexTPC(0.),
56   fMaxDCAXY(0.), fMaxDCAXYTPC(0.),
57   fMaxDCAZ(0.), fMaxDCAZTPC(0.),
58   fMaxDCA3D(0.), fMaxDCA3DTPC(0.),
59   fMaxConstrainChi2(0.), fMinTPCdEdxPoints(0),
60   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
61   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
62   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
63   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
64   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
65   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
66   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
67   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
68   fMaxConstrainChi2Flag(kFALSE),
69   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
70   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
71   fPointOnSPDLayersFlag(0), fPointOnSDDLayersFlag(0), fPointOnSSDLayersFlag(0),
72   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
73   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
74   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
75   fMinTPCdEdxPointsFlag(kFALSE),
76   fPtDependentDcaXY(0), fPtDependentDcaXYFlag(kFALSE), fNSigmaDCAXY(0),
77   fFunctionProbabilityFlag(kFALSE), 
78   fNSigma(0), fNRatio(0.),
79   fElectronFunction(0), fMuonFunction(0),
80   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
81   fDebugMode(kFALSE), fListVertexQA(0) {
82   //Default constructor
83   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
84   /*for(Int_t i = 0; i < 24; i++) {
85     fdEdxMean[i] = 0.0;
86     fdEdxSigma[i] = 0.0;
87     }*/
88
89   fListVertexQA = new TList();
90   fListVertexQA->SetName("fListVertexQA");
91   TH1F *gHistVx = new TH1F("gHistVx",
92                            "Vx distribution;V_{x} [cm];Entries",
93                            100,-5.,5.);
94   gHistVx->SetFillColor(kRed-2);
95   fListVertexQA->Add(gHistVx);
96   TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
97                                    "Vx distribution;V_{x} [cm];Entries",
98                                    100,-5.,5.);
99   fListVertexQA->Add(gHistVxAccepted);
100   TH1F *gHistVy = new TH1F("gHistVy",
101                            "Vy distribution;V_{y} [cm];Entries",
102                            100,-5.,5.);
103   gHistVy->SetFillColor(kRed-2);
104   fListVertexQA->Add(gHistVy);
105   TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
106                                    "Vy distribution;V_{y} [cm];Entries",
107                                    100,-5.,5.);
108   fListVertexQA->Add(gHistVyAccepted);
109   TH1F *gHistVz = new TH1F("gHistVz",
110                            "Vz distribution;V_{z} [cm];Entries",
111                            100,-25.,25.);
112   gHistVz->SetFillColor(kRed-2);
113   fListVertexQA->Add(gHistVz);
114   TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
115                                    "Vz distribution;V_{z} [cm];Entries",
116                                    100,-25.,25.);
117   fListVertexQA->Add(gHistVzAccepted);
118
119   TH1F *gHistNumberOfContributors = new TH1F("gHistNumberOfContributors",
120                                              "Number of contributors;N_{contr.};Entries",
121                                              100,0.,100.);
122   fListVertexQA->Add(gHistNumberOfContributors);
123
124
125 }
126
127 //____________________________________________________________________//
128 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
129   //Default destructor
130   if(fElectronFunction) delete fElectronFunction;
131   if(fMuonFunction) delete fMuonFunction;
132   if(fPionFunction) delete fPionFunction;
133   if(fKaonFunction) delete fKaonFunction;
134   if(fProtonFunction) delete fProtonFunction;
135   if(fListVertexQA) delete fListVertexQA;
136   if(fPtDependentDcaXY) delete fPtDependentDcaXY;
137   if(fPhysicsSelection) delete fPhysicsSelection;
138 }
139
140 //____________________________________________________________________//
141 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
142   //Return the a priori probs
143   Double_t partFrac=0;
144   if(fFunctionProbabilityFlag) {
145     if(i == 0) partFrac = fElectronFunction->Eval(p);
146     if(i == 1) partFrac = fMuonFunction->Eval(p);
147     if(i == 2) partFrac = fPionFunction->Eval(p);
148     if(i == 3) partFrac = fKaonFunction->Eval(p);
149     if(i == 4) partFrac = fProtonFunction->Eval(p);
150   }
151   else partFrac = fPartFrac[i];
152
153   return partFrac;
154 }
155
156 //____________________________________________________________________//
157 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
158   // Checks if the track is outside the analyzed y-Pt phase space
159   Double_t gP = 0.0, gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
160   Double_t eta = 0.0;
161
162   if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid) || (fProtonAnalysisMode == kFullHybrid)) {
163     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
164     if(!tpcTrack) {
165       gP = 0.0; gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
166     }
167     else {
168       gP = tpcTrack->P();
169       gPt = tpcTrack->Pt();
170       gPx = tpcTrack->Px();
171       gPy = tpcTrack->Py();
172       gPz = tpcTrack->Pz();
173       eta = tpcTrack->Eta();
174     }
175   }//standalone TPC or Hybrid TPC approaches
176   else {
177     gP = track->P();
178     gPt = track->Pt();
179     gPx = track->Px();
180     gPy = track->Py();
181     gPz = track->Pz();
182     eta = track->Eta();
183   }
184   
185   if((gPt < fMinY) || (gPt > fMaxY)) {
186       if(fDebugMode)
187         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
188       return kFALSE;
189   }
190   if((gP < fMinY) || (gP > fMaxY)) {
191       if(fDebugMode)
192         Printf("IsInPhaseSpace: Track rejected because it has a P value of %lf (accepted interval: %lf - %lf)",gP,fMinY,fMaxY);
193       return kFALSE;
194   }
195   if(fAnalysisEtaMode) {
196     if((eta < fMinX) || (eta > fMaxX)) {
197       if(fDebugMode)
198         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
199       return kFALSE;
200     }
201   }
202   else {
203     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
204       if(fDebugMode)
205         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
206       return kFALSE;
207     }
208   }
209
210   return kTRUE;
211 }
212
213 //____________________________________________________________________//
214 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDtrack* track) {
215   // Checks if the track is excluded from the cuts
216   //Int_t  fIdxInt[200];
217   //Int_t nClustersITS = track->GetITSclusters(fIdxInt);
218   //Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
219   Int_t nClustersITS = track->GetITSclusters(0x0);
220   Int_t nClustersTPC = track->GetTPCclusters(0x0);
221
222   Float_t chi2PerClusterITS = -1;
223   if (nClustersITS!=0)
224     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
225   Float_t chi2PerClusterTPC = -1;
226   if (nClustersTPC!=0)
227     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
228
229   Double_t extCov[15];
230   track->GetExternalCovariance(extCov);
231
232   if(fPointOnSPDLayersFlag) {
233     if((!track->HasPointOnITSLayer(0))&&(!track->HasPointOnITSLayer(1))) {
234       if(fDebugMode)
235         Printf("IsAccepted: Track rejected because it doesn't have a point on either SPD layers");
236       return kFALSE;
237     }
238   }
239   if(fPointOnSDDLayersFlag) {
240     if((!track->HasPointOnITSLayer(2))&&(!track->HasPointOnITSLayer(3))) {
241       if(fDebugMode)
242         Printf("IsAccepted: Track rejected because it doesn't have a point on either SDD layers");
243       return kFALSE;
244     }
245   }
246   if(fPointOnSSDLayersFlag) {
247     if((!track->HasPointOnITSLayer(4))&&(!track->HasPointOnITSLayer(5))) {
248       if(fDebugMode)
249         Printf("IsAccepted: Track rejected because it doesn't have a point on either SSD layers");
250       return kFALSE;
251     }
252   }
253   if(fPointOnITSLayer1Flag) {
254     if(!track->HasPointOnITSLayer(0)) {
255       if(fDebugMode)
256         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
257       return kFALSE;
258     }
259   }
260   if(fPointOnITSLayer2Flag) {
261     if(!track->HasPointOnITSLayer(1)) {
262       if(fDebugMode)
263         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
264       return kFALSE;
265     }
266   }
267   if(fPointOnITSLayer3Flag) {
268     if(!track->HasPointOnITSLayer(2)) {
269       if(fDebugMode)
270         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
271       return kFALSE;
272     }
273   }
274   if(fPointOnITSLayer4Flag) {
275     if(!track->HasPointOnITSLayer(3)) {
276       if(fDebugMode)
277         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
278       return kFALSE;
279     }
280   }
281   if(fPointOnITSLayer5Flag) {
282     if(!track->HasPointOnITSLayer(4)) {
283       if(fDebugMode)
284         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
285       return kFALSE;
286     }
287   }
288   if(fPointOnITSLayer6Flag) {
289     if(!track->HasPointOnITSLayer(5)) {
290       if(fDebugMode)
291         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
292       return kFALSE;
293     }
294   }
295   if(fMinITSClustersFlag) {
296     if(nClustersITS < fMinITSClusters) {
297       if(fDebugMode)
298         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
299       return kFALSE;
300     }
301   }
302   if(fMaxChi2PerITSClusterFlag) {
303     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
304       if(fDebugMode)
305         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
306       return kFALSE; 
307     }
308   }
309   if(fMinTPCClustersFlag) {
310     if(nClustersTPC < fMinTPCClusters) {
311       if(fDebugMode)
312         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
313       return kFALSE;
314     }
315   }
316   if(fMaxChi2PerTPCClusterFlag) {
317     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
318       if(fDebugMode)
319         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
320       return kFALSE; 
321     }
322   }
323   if(fMaxCov11Flag) {
324     if(extCov[0] > fMaxCov11) {
325       if(fDebugMode)
326         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
327       return kFALSE;
328     }
329   }
330   if(fMaxCov22Flag) {
331     if(extCov[2] > fMaxCov22) {
332       if(fDebugMode)
333         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
334       return kFALSE;
335     }
336   }
337   if(fMaxCov33Flag) {
338     if(extCov[5] > fMaxCov33) {
339       if(fDebugMode)
340         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
341       return kFALSE;
342     }
343   }
344   if(fMaxCov44Flag) {
345     if(extCov[9] > fMaxCov44) {
346       if(fDebugMode)
347         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
348       return kFALSE;
349     }
350   }
351   if(fMaxCov55Flag) {
352     if(extCov[14] > fMaxCov55) {
353       if(fDebugMode)
354         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
355       return kFALSE;
356     }
357   }
358   if(fMinTPCdEdxPointsFlag) {
359     if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
360       if(fDebugMode)
361         Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
362       return kFALSE;
363     }
364   }
365   if(fITSRefitFlag) {
366     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
367       if(fDebugMode)
368         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
369       return kFALSE;
370     }
371   }
372   if(fTPCRefitFlag) {
373     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
374       if(fDebugMode)
375         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
376       return kFALSE;
377     }
378   }
379   if(fESDpidFlag) {
380     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
381       if(fDebugMode)
382         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
383       return kFALSE;
384     }
385   }
386   if(fTPCpidFlag) {
387     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
388       if(fDebugMode)
389         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
390       return kFALSE;
391     }
392   }
393   if(fTOFpidFlag) {
394     if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
395       if(fDebugMode)
396         Printf("IsAccepted: Track rejected because it has no TOF pid flag");
397       return kFALSE;
398     }
399   }
400
401   return kTRUE;
402 }
403
404 //____________________________________________________________________//
405 void AliProtonAnalysisBase::SetPtDependentDCAxy(Int_t nSigma, 
406                                                 Double_t p0, 
407                                                 Double_t p1, 
408                                                 Double_t p2) {
409   //Pt dependent dca cut in xy
410   fNSigmaDCAXY = nSigma;
411   fPtDependentDcaXY = new TF1("fPtDependentDcaXY","[0]+[1]/x^[2]",0.1,10.1);
412   fPtDependentDcaXY->SetParameter(0,p0); 
413   fPtDependentDcaXY->SetParameter(1,p1);
414   fPtDependentDcaXY->SetParameter(2,p2);
415   fPtDependentDcaXYFlag = kTRUE;
416 }
417
418 //____________________________________________________________________//
419 Bool_t AliProtonAnalysisBase::IsPrimary(AliESDEvent *esd,
420                                         const AliESDVertex *vertex, 
421                                         AliESDtrack* track) {
422   // Checks if the track is a primary-like candidate
423   const Double_t kMicrometer2Centimeter = 0.0001;
424   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
425   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
426   Double_t dca3D = 0.0;
427   Float_t dcaXY = 0.0, dcaZ = 0.0;
428
429   if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
430     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
431     if(!tpcTrack) {
432       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
433       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
434       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
435     }
436     else {
437       gPt = tpcTrack->Pt();
438       gPx = tpcTrack->Px();
439       gPy = tpcTrack->Py();
440       gPz = tpcTrack->Pz();
441       tpcTrack->PropagateToDCA(vertex,
442                                esd->GetMagneticField(),
443                                100.,dca,cov);
444     }
445   }//standalone TPC or hybrid TPC approaches
446   if(fProtonAnalysisMode == kFullHybrid) {
447     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
448     AliExternalTrackParam cParam;
449     if(!tpcTrack) {
450       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
451       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
452       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
453     }
454     else {
455       gPt = tpcTrack->Pt();
456       gPx = tpcTrack->Px();
457       gPy = tpcTrack->Py();
458       gPz = tpcTrack->Pz();
459       track->RelateToVertex(vertex,
460                             esd->GetMagneticField(),
461                             100.,&cParam);
462       track->GetImpactParameters(dcaXY,dcaZ);
463       dca[0] = dcaXY; dca[1] = dcaZ;
464     }
465   }//standalone TPC or hybrid TPC approaches
466   else {
467     gPt = track->Pt();
468     gPx = track->Px();
469     gPy = track->Py();
470     gPz = track->Pz();
471     track->PropagateToDCA(vertex,
472                           esd->GetMagneticField(),
473                           100.,dca,cov);
474   }
475   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
476                       TMath::Power(dca[1],2));
477      
478   if(fMaxSigmaToVertexFlag) {
479     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
480       if(fDebugMode)
481         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
482       return kFALSE;
483     }
484   }
485   if(fMaxSigmaToVertexTPCFlag) {
486     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
487       if(fDebugMode)
488         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
489       return kFALSE;
490     }
491   }
492   if(fMaxDCAXYFlag) { 
493     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
494       if(fDebugMode)
495         Printf("IsPrimary: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
496       return kFALSE;
497     }
498   }
499   if(fMaxDCAXYTPCFlag) { 
500     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
501       if(fDebugMode)
502         Printf("IsPrimary: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
503       return kFALSE;
504     }
505   }
506   if(fMaxDCAZFlag) { 
507     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
508       if(fDebugMode)
509         Printf("IsPrimary: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
510       return kFALSE;
511     }
512   }
513   if(fMaxDCAZTPCFlag) { 
514     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
515       if(fDebugMode)
516         Printf("IsPrimary: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
517       return kFALSE;
518     }
519   }
520   if(fMaxDCA3DFlag) { 
521     if(TMath::Abs(dca3D) > fMaxDCA3D) {
522       if(fDebugMode)
523         Printf("IsPrimary: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
524       return kFALSE;
525     }
526   }
527   if(fMaxDCA3DTPCFlag) { 
528     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
529       if(fDebugMode)
530         Printf("IsPrimary: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
531       return kFALSE;
532     }
533   }
534   if(fMaxConstrainChi2Flag) {
535     if(track->GetConstrainedChi2() > 0) 
536       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
537       if(fDebugMode)
538         Printf("IsPrimary: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
539       return kFALSE;
540       }
541   }
542   if(fPtDependentDcaXYFlag) {
543     if(TMath::Abs(dca[0]) > kMicrometer2Centimeter*fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt)) {
544       if(fDebugMode)
545         Printf("IsPrimary: Track rejected because it has a value of the dca(xy) higher than the %d sigma pt dependent cut: %lf (max. requested: %lf)",fNSigmaDCAXY,TMath::Abs(dca[0]),fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt));
546       return kFALSE;
547     }
548   }
549
550   return kTRUE;
551 }
552
553 //____________________________________________________________________//
554 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
555   // Calculates the number of sigma to the vertex.
556   Float_t b[2];
557   Float_t bRes[2];
558   Float_t bCov[3];
559   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid)&&(fProtonAnalysisMode != kFullHybrid))
560     esdTrack->GetImpactParametersTPC(b,bCov);
561   else
562     esdTrack->GetImpactParameters(b,bCov);
563   
564   if (bCov[0]<=0 || bCov[2]<=0) {
565     //AliDebug(1, "Estimated b resolution lower or equal zero!");
566     bCov[0]=0; bCov[2]=0;
567   }
568   bRes[0] = TMath::Sqrt(bCov[0]);
569   bRes[1] = TMath::Sqrt(bCov[2]);
570   
571   if (bRes[0] == 0 || bRes[1] ==0) return -1;
572   
573   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
574   
575   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
576   
577   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
578   
579   return d;
580 }
581
582 //____________________________________________________________________//
583 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx, 
584                                          Double_t gPy, 
585                                          Double_t gPz) const {
586   //returns the rapidity of the proton - to be removed
587   Double_t fMass = 9.38270000000000048e-01;
588   
589   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) + 
590                            TMath::Power(gPy,2) + 
591                            TMath::Power(gPz,2));
592   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
593   Double_t y = -999;
594   if(energy != gPz) 
595     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
596
597   return y;
598 }
599
600 //________________________________________________________________________
601 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd, 
602                                                      AnalysisMode mode,
603                                                      Double_t gVxMax,
604                                                      Double_t gVyMax,
605                                                      Double_t gVzMax) {
606   // Get the vertex from the ESD and returns it if the vertex is valid
607   // Second argument decides which vertex is used (this selects
608   // also the quality criteria that are applied)
609   const AliESDVertex* vertex = 0;
610   if((mode == kHybrid)||(mode == kFullHybrid))
611     vertex = esd->GetPrimaryVertexSPD();
612   else if(mode == kTPC){
613     Double_t kBz = esd->GetMagneticField();
614     AliVertexerTracks vertexer(kBz);
615     vertexer.SetTPCMode();
616     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
617     esd->SetPrimaryVertexTPC(vTPC);
618     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
619       AliESDtrack *t = esd->GetTrack(i);
620       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
621     }
622     delete vTPC;
623     vertex = esd->GetPrimaryVertexTPC();
624   }
625   else if(mode == kGlobal)
626     vertex = esd->GetPrimaryVertex();
627   else
628     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
629   
630   if(!vertex) {
631     if(fDebugMode)
632       Printf("GetVertex: Event rejected because there is no valid vertex object");
633     return 0;
634   }
635   
636   // check Ncontributors
637   if(vertex->GetNContributors() <= 0) {
638     if(fDebugMode)
639       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
640     return 0;
641   }
642   
643   // check resolution
644   Double_t zRes = vertex->GetZRes();
645   if(zRes == 0) {
646     if(fDebugMode)
647       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
648     return 0;
649   }
650   ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
651   ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
652   ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
653
654   //check position
655   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
656     if(fDebugMode)
657       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
658     return 0;
659   }
660   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
661     if(fDebugMode)
662       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
663     return 0;
664   }
665   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
666     if(fDebugMode)
667       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
668     return 0;
669   }
670   ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
671   ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
672   ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
673   ((TH1F *)(fListVertexQA->At(6)))->Fill(vertex->GetNContributors());
674
675   //check number of contributors
676   if(fMinNumOfContributors > 0) {
677     if(fMinNumOfContributors > vertex->GetNContributors()) {
678       if(fDebugMode)
679         Printf("GetVertex: Event rejected because it has %d number of contributors (requested minimum: %d)",vertex->GetNContributors(),fMinNumOfContributors);
680       
681       return 0;
682     }
683   }
684   
685   return vertex;
686 }
687
688 //________________________________________________________________________
689 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd, 
690                                                TriggerMode trigger) {
691   // check if the event was triggered
692   ULong64_t triggerMask = esd->GetTriggerMask();
693   TString firedTriggerClass = esd->GetFiredTriggerClasses();
694
695   if(fAnalysisMC) {
696     // definitions from p-p.cfg
697     ULong64_t spdFO = (1 << 14);
698     ULong64_t v0left = (1 << 11);
699     ULong64_t v0right = (1 << 12);
700     
701     switch (trigger) {
702     case kMB1: {
703       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
704         return kTRUE;
705       break;
706     }
707     case kMB2: {
708       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
709         return kTRUE;
710       break;
711     }
712     case kSPDFASTOR: {
713       if (triggerMask & spdFO)
714         return kTRUE;
715       break;
716     }
717     }//switch
718   }
719   else {
720     if(kUseOnlineTrigger) {
721       if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
722         return kTRUE;
723     }
724     else if(!kUseOnlineTrigger)
725       return kTRUE;
726   }
727
728   return kFALSE;
729 }
730
731 //________________________________________________________________________
732 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
733   // return the list of cuts and their cut values for reference
734   TLatex l;
735   l.SetTextAlign(12);
736   l.SetTextSize(0.04);
737
738   TString listOfCuts;
739
740   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
741   c->SetFillColor(10); c->SetHighLightColor(41);
742   c->Divide(3,2);
743
744   c->cd(1)->SetFillColor(10);
745   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
746
747   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
748   l.DrawLatex(0.1,0.82,listOfCuts.Data());
749   listOfCuts = "Analysis mode: ";
750   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone"; 
751   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC"; 
752   if(fProtonAnalysisMode == kFullHybrid) listOfCuts += "Full Hybrid TPC"; 
753   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking"; 
754   l.DrawLatex(0.1,0.74,listOfCuts.Data());
755   listOfCuts = "Trigger mode: "; 
756   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1"; 
757   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2"; 
758   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)"; 
759   l.DrawLatex(0.1,0.66,listOfCuts.Data());
760   listOfCuts = "PID mode: "; 
761   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
762   if(fProtonPIDMode == kRatio) {
763     listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.}) > ";
764     listOfCuts += fNRatio;
765   }
766   if(fProtonPIDMode == kSigma) {
767     listOfCuts += "N_{#sigma} area: "; listOfCuts += fNSigma;
768     listOfCuts += " #sigma";
769   }
770   //if(fProtonPIDMode == kSigma2) {
771   //listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
772   //listOfCuts += " #sigma";
773   //}
774   l.DrawLatex(0.1,0.58,listOfCuts.Data());
775   listOfCuts = "Accepted vertex diamond: "; 
776   l.DrawLatex(0.1,0.52,listOfCuts.Data());
777   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
778   l.DrawLatex(0.6,0.52,listOfCuts.Data());
779   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
780   l.DrawLatex(0.6,0.45,listOfCuts.Data());
781   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
782   l.DrawLatex(0.6,0.38,listOfCuts.Data());
783   listOfCuts = "N_{contributors} > "; listOfCuts += fMinNumOfContributors; 
784   l.DrawLatex(0.6,0.31,listOfCuts.Data());
785   listOfCuts = "Phase space: "; 
786   l.DrawLatex(0.1,0.2,listOfCuts.Data());
787   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";  
788   else listOfCuts = "|y| < ";
789   listOfCuts += TMath::Abs(fMinX); 
790   l.DrawLatex(0.35,0.2,listOfCuts.Data());
791   listOfCuts = "N_{bins} = ";
792   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
793   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
794   l.DrawLatex(0.35,0.15,listOfCuts.Data());
795   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
796   listOfCuts += fMaxY; listOfCuts += "GeV/c"; 
797   l.DrawLatex(0.35,0.1,listOfCuts.Data());
798   listOfCuts = "N_{bins} = ";
799   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
800   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
801   l.DrawLatex(0.35,0.05,listOfCuts.Data());
802
803   c->cd(2)->SetFillColor(2);
804   l.DrawLatex(0.3,0.95,"ITS related cuts");
805   listOfCuts = "Request a cluster on either of the SPD layers: "; 
806   listOfCuts += fPointOnSPDLayersFlag;
807   l.DrawLatex(0.1,0.9,listOfCuts.Data());
808   listOfCuts = "Request a cluster on either of the SDD layers: "; 
809   listOfCuts += fPointOnSDDLayersFlag;
810   l.DrawLatex(0.1,0.83,listOfCuts.Data());
811   listOfCuts = "Request a cluster on either of the SSD layers: "; 
812   listOfCuts += fPointOnSSDLayersFlag;
813   l.DrawLatex(0.1,0.76,listOfCuts.Data());
814
815   listOfCuts = "Request a cluster on SPD1: "; 
816   listOfCuts += fPointOnITSLayer1Flag;
817   l.DrawLatex(0.1,0.69,listOfCuts.Data());
818   listOfCuts = "Request a cluster on SPD2: "; 
819   listOfCuts += fPointOnITSLayer2Flag;
820   l.DrawLatex(0.1,0.62,listOfCuts.Data());
821   listOfCuts = "Request a cluster on SDD1: "; 
822   listOfCuts += fPointOnITSLayer3Flag;
823   l.DrawLatex(0.1,0.55,listOfCuts.Data());
824   listOfCuts = "Request a cluster on SDD2: "; 
825   listOfCuts += fPointOnITSLayer4Flag;
826   l.DrawLatex(0.1,0.48,listOfCuts.Data());
827   listOfCuts = "Request a cluster on SSD1: "; 
828   listOfCuts += fPointOnITSLayer5Flag;
829   l.DrawLatex(0.1,0.41,listOfCuts.Data());
830   listOfCuts = "Request a cluster on SSD2: "; 
831   listOfCuts += fPointOnITSLayer6Flag; 
832   l.DrawLatex(0.1,0.34,listOfCuts.Data());  
833   listOfCuts = "Minimum number of ITS clusters: ";
834   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
835   else listOfCuts += "Not used";
836   l.DrawLatex(0.1,0.27,listOfCuts.Data());
837   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
838   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster; 
839   else listOfCuts += "Not used";
840   l.DrawLatex(0.1,0.2,listOfCuts.Data());
841
842   c->cd(3)->SetFillColor(3);
843   l.DrawLatex(0.3,0.9,"TPC related cuts");
844   listOfCuts = "Minimum number of TPC clusters: ";
845   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
846   else listOfCuts += "Not used";
847   l.DrawLatex(0.1,0.7,listOfCuts.Data());
848   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
849   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
850   else listOfCuts += "Not used";
851   l.DrawLatex(0.1,0.5,listOfCuts.Data());
852   listOfCuts = "Minimum number of TPC points for the dE/dx: ";
853   if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
854   else listOfCuts += "Not used";
855   l.DrawLatex(0.1,0.3,listOfCuts.Data());
856
857   c->cd(4)->SetFillColor(4);
858   l.DrawLatex(0.3,0.9,"Tracking related cuts");
859   listOfCuts = "Maximum cov11: ";
860   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
861   else listOfCuts += "Not used";
862   l.DrawLatex(0.1,0.75,listOfCuts.Data());
863   listOfCuts = "Maximum cov22: ";
864   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
865   else listOfCuts += "Not used";
866   l.DrawLatex(0.1,0.6,listOfCuts.Data());
867   listOfCuts = "Maximum cov33: ";
868   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
869   else listOfCuts += "Not used";
870   l.DrawLatex(0.1,0.45,listOfCuts.Data());
871   listOfCuts = "Maximum cov44: ";
872   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
873   else listOfCuts += "Not used";
874   l.DrawLatex(0.1,0.3,listOfCuts.Data());
875   listOfCuts = "Maximum cov55: ";
876   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
877   else listOfCuts += "Not used";
878   l.DrawLatex(0.1,0.15,listOfCuts.Data());
879
880   c->cd(5)->SetFillColor(5);
881   l.DrawLatex(0.3,0.9,"DCA related cuts");
882   listOfCuts = "Maximum sigma to vertex: ";
883   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
884   else listOfCuts += "Not used";
885   l.DrawLatex(0.1,0.8,listOfCuts.Data());
886   listOfCuts = "Maximum sigma to vertex (TPC): ";
887   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
888   else listOfCuts += "Not used";
889   l.DrawLatex(0.1,0.72,listOfCuts.Data());
890   listOfCuts = "Maximum DCA in xy: ";
891   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
892   else listOfCuts += "Not used";
893   l.DrawLatex(0.1,0.64,listOfCuts.Data());
894   listOfCuts = "Maximum DCA in z: ";
895   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
896   else listOfCuts += "Not used";
897   l.DrawLatex(0.1,0.56,listOfCuts.Data());
898   listOfCuts = "Maximum DCA in 3D: ";
899   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
900   else listOfCuts += "Not used";
901   l.DrawLatex(0.1,0.48,listOfCuts.Data());
902   listOfCuts = "Maximum DCA in xy (TPC): ";
903   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
904   else listOfCuts += "Not used";
905   l.DrawLatex(0.1,0.4,listOfCuts.Data());
906   listOfCuts = "Maximum DCA in z (TPC): ";
907   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
908   else listOfCuts += "Not used";
909   l.DrawLatex(0.1,0.32,listOfCuts.Data());
910   listOfCuts = "Maximum DCA in 3D (TPC): ";
911   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
912   else listOfCuts += "Not used";
913   l.DrawLatex(0.1,0.24,listOfCuts.Data());
914   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
915   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
916   else listOfCuts += "Not used";
917   l.DrawLatex(0.1,0.16,listOfCuts.Data());
918
919   c->cd(6)->SetFillColor(6);
920   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
921   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
922   l.DrawLatex(0.1,0.7,listOfCuts.Data());
923   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag; 
924   l.DrawLatex(0.1,0.5,listOfCuts.Data());
925   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
926   l.DrawLatex(0.1,0.3,listOfCuts.Data());
927   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
928   l.DrawLatex(0.1,0.1,listOfCuts.Data());
929
930   return c;
931 }
932
933 //________________________________________________________________________
934 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
935   //Function that checks if a track is a proton
936   Double_t probability[5] = {0.,0.,0.,0.,0.};
937   Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
938   Long64_t fParticleType = 0;
939  
940   //Bayesian approach for the PID
941   if(fProtonPIDMode == kBayesian) {
942     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)||(fProtonAnalysisMode == kFullHybrid)) {
943       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
944       if(tpcTrack) {
945         gPt = tpcTrack->Pt();
946         gP = tpcTrack->P();
947         track->GetTPCpid(probability);
948       }
949       else return 0;
950     }//TPC standalone or Hybrid TPC
951     else if(fProtonAnalysisMode == kGlobal) {
952       gPt = track->Pt();
953       gP = track->P();
954       track->GetESDpid(probability);
955     }//Global tracking    
956     
957     Double_t rcc = 0.0;
958     for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
959       rcc += probability[i]*GetParticleFraction(i,gP);
960     if(rcc != 0.0) {
961       Double_t w[5];
962       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
963         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
964       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
965     }
966     if(fParticleType == 4)
967       return kTRUE;
968   }//Bayesian pid
969   //Ratio of the measured over the theoretical dE/dx a la STAR
970   else if(fProtonPIDMode == kRatio) {
971     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
972     if(tpcTrack) {
973       gPt = tpcTrack->Pt();
974       gP = track->GetInnerParam()->P();
975       gEta = tpcTrack->Eta();
976     }
977     Double_t fAlephParameters[5] = {0.,0.,0.,0.,0.};
978     if(fAnalysisMC) {
979       fAlephParameters[0] = 2.15898e+00/50.;
980       fAlephParameters[1] = 1.75295e+01;
981       fAlephParameters[2] = 3.40030e-09;
982       fAlephParameters[3] = 1.96178e+00;
983       fAlephParameters[4] = 3.91720e+00;
984     }
985     else {
986       fAlephParameters[0] = 0.0283086;
987       fAlephParameters[1] = 2.63394e+01;
988       fAlephParameters[2] = 5.04114e-11;
989       fAlephParameters[3] = 2.12543e+00;
990       fAlephParameters[4] = 4.88663e+00;
991     }
992     
993     AliTPCPIDResponse *tpcResponse = new AliTPCPIDResponse();
994     tpcResponse->SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
995
996     Double_t normalizeddEdx = -10.;
997     if((track->GetTPCsignal() > 0.0) && (tpcResponse->GetExpectedSignal(gP,AliPID::kProton) > 0.0))
998       normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse->GetExpectedSignal(gP,AliPID::kProton));
999
1000     delete tpcResponse;
1001
1002     if(normalizeddEdx >= fNRatio)
1003       return kTRUE;
1004   }//kRatio PID mode
1005   //Definition of an N-sigma area around the dE/dx vs P band
1006   else if(fProtonPIDMode == kSigma) {
1007     Double_t fAlephParameters[5];
1008     if(fAnalysisMC) {
1009       fAlephParameters[0] = 2.15898e+00/50.;
1010       fAlephParameters[1] = 1.75295e+01;
1011       fAlephParameters[2] = 3.40030e-09;
1012       fAlephParameters[3] = 1.96178e+00;
1013       fAlephParameters[4] = 3.91720e+00;
1014     }
1015     else {
1016       fAlephParameters[0] = 0.0283086;
1017       fAlephParameters[1] = 2.63394e+01;
1018       fAlephParameters[2] = 5.04114e-11;
1019       fAlephParameters[3] = 2.12543e+00;
1020       fAlephParameters[4] = 4.88663e+00;
1021     }
1022
1023     Double_t nsigma = 100.0;
1024     AliTPCPIDResponse *tpcResponse = new AliTPCPIDResponse();
1025     tpcResponse->SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
1026     
1027     Double_t mom = track->GetP();
1028     const AliExternalTrackParam *in = track->GetInnerParam();
1029     if (in)
1030       mom = in->GetP();
1031
1032     nsigma = TMath::Abs(tpcResponse->GetNumberOfSigmas(mom,track->GetTPCsignal(),track->GetTPCsignalN(),AliPID::kProton));
1033   
1034     delete tpcResponse;
1035     if(nsigma <= fNSigma) 
1036       return kTRUE;
1037   }//kSigma PID method
1038
1039   return kFALSE;
1040 }
1041
1042
1043