Adding READMEeve.txt page
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / 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
26 #include <AliExternalTrackParam.h>
27 #include <AliESDEvent.h>
28 #include <AliPID.h>
29 #include <AliVertexerTracks.h>
30
31 class AliLog;
32 class AliESDVertex;
33
34 #include "AliProtonAnalysisBase.h"
35
36 ClassImp(AliProtonAnalysisBase)
37
38 //____________________________________________________________________//
39 AliProtonAnalysisBase::AliProtonAnalysisBase() : 
40   TObject(),  fProtonAnalysisLevel("ESD"),
41   fTriggerMode(kMB2), fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
42   fAnalysisEtaMode(kFALSE),
43   fVxMax(100.), fVyMax(100.), fVzMax(100.),
44   fNBinsX(0), fMinX(0), fMaxX(0),
45   fNBinsY(0), fMinY(0), fMaxY(0),
46   fMinTPCClusters(0), fMinITSClusters(0),
47   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
48   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
49   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
50   fMaxDCAXY(0), fMaxDCAXYTPC(0),
51   fMaxDCAZ(0), fMaxDCAZTPC(0),
52   fMaxDCA3D(0), fMaxDCA3DTPC(0),
53   fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
54   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
55   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
56   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
57   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
58   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
59   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
60   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
61   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
62   fMaxConstrainChi2Flag(kFALSE),
63   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
64   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
65   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
66   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
67   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
68   fMinTPCdEdxPointsFlag(kFALSE),
69   fFunctionProbabilityFlag(kFALSE), 
70   fNSigma(0),
71   fElectronFunction(0), fMuonFunction(0),
72   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
73   fDebugMode(kFALSE) {
74   //Default constructor
75   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
76   for(Int_t i = 0; i < 24; i++) {
77     fdEdxMean[i] = 0.0;
78     fdEdxSigma[i] = 0.0;
79   }
80 }
81
82 //____________________________________________________________________//
83 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
84   //Default destructor
85   if(fElectronFunction) delete fElectronFunction;
86   if(fMuonFunction) delete fMuonFunction;
87   if(fPionFunction) delete fPionFunction;
88   if(fKaonFunction) delete fKaonFunction;
89   if(fProtonFunction) delete fProtonFunction;
90 }
91
92 //____________________________________________________________________//
93 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
94   //Return the a priori probs
95   Double_t partFrac=0;
96   if(fFunctionProbabilityFlag) {
97     if(i == 0) partFrac = fElectronFunction->Eval(p);
98     if(i == 1) partFrac = fMuonFunction->Eval(p);
99     if(i == 2) partFrac = fPionFunction->Eval(p);
100     if(i == 3) partFrac = fKaonFunction->Eval(p);
101     if(i == 4) partFrac = fProtonFunction->Eval(p);
102   }
103   else partFrac = fPartFrac[i];
104
105   return partFrac;
106 }
107
108 //____________________________________________________________________//
109 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
110   // Checks if the track is outside the analyzed y-Pt phase space
111   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
112   Double_t eta = 0.0;
113
114   if(fProtonAnalysisMode == kTPC) {
115     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
116     if(!tpcTrack) {
117       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
118     }
119     else {
120       gPt = tpcTrack->Pt();
121       gPx = tpcTrack->Px();
122       gPy = tpcTrack->Py();
123       gPz = tpcTrack->Pz();
124       eta = tpcTrack->Eta();
125     }
126   }
127   else {
128     gPt = track->Pt();
129     gPx = track->Px();
130     gPy = track->Py();
131     gPz = track->Pz();
132     eta = track->Eta();
133   }
134   
135   if((gPt < fMinY) || (gPt > fMaxY)) {
136       if(fDebugMode)
137         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
138       return kFALSE;
139   }
140   if(fAnalysisEtaMode) {
141     if((eta < fMinX) || (eta > fMaxX)) {
142       if(fDebugMode)
143         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
144       return kFALSE;
145     }
146   }
147   else {
148     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
149       if(fDebugMode)
150         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
151       return kFALSE;
152     }
153   }
154
155   return kTRUE;
156 }
157
158 //____________________________________________________________________//
159 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
160                                          const AliESDVertex *vertex, 
161                                          AliESDtrack* track) {
162   // Checks if the track is excluded from the cuts
163   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
164   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
165   Double_t dca3D = 0.0;
166   
167   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid)) {
168     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
169     if(!tpcTrack) {
170       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
171       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
172       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
173     }
174     else {
175       gPt = tpcTrack->Pt();
176       gPx = tpcTrack->Px();
177       gPy = tpcTrack->Py();
178       gPz = tpcTrack->Pz();
179       tpcTrack->PropagateToDCA(vertex,
180                                esd->GetMagneticField(),
181                                100.,dca,cov);
182     }
183   }
184   else if(fProtonAnalysisMode == kHybrid) {
185      AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
186     if(!tpcTrack) {
187       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
188       dca[0] = -100.; dca[1] = -100.;
189       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
190     }
191     else {
192       gPt = tpcTrack->Pt();
193       gPx = tpcTrack->Px();
194       gPy = tpcTrack->Py();
195       gPz = tpcTrack->Pz();
196       tpcTrack->PropagateToDCA(vertex,
197                                esd->GetMagneticField(),
198                                100.,dca,cov);
199     }
200   }
201   else{
202     gPt = track->Pt();
203     gPx = track->Px();
204     gPy = track->Py();
205     gPz = track->Pz();
206     track->PropagateToDCA(vertex,
207                           esd->GetMagneticField(),
208                           100.,dca,cov);
209   }
210   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
211                       TMath::Power(dca[1],2));
212      
213   Int_t  fIdxInt[200];
214   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
215   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
216
217   Float_t chi2PerClusterITS = -1;
218   if (nClustersITS!=0)
219     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
220   Float_t chi2PerClusterTPC = -1;
221   if (nClustersTPC!=0)
222     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
223
224   Double_t extCov[15];
225   track->GetExternalCovariance(extCov);
226
227   if(fPointOnITSLayer1Flag) {
228     if(!track->HasPointOnITSLayer(0)) {
229       if(fDebugMode)
230         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
231       return kFALSE;
232     }
233   }
234   if(fPointOnITSLayer2Flag) {
235     if(!track->HasPointOnITSLayer(1)) {
236       if(fDebugMode)
237         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
238       return kFALSE;
239     }
240   }
241   if(fPointOnITSLayer3Flag) {
242     if(!track->HasPointOnITSLayer(2)) {
243       if(fDebugMode)
244         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
245       return kFALSE;
246     }
247   }
248   if(fPointOnITSLayer4Flag) {
249     if(!track->HasPointOnITSLayer(3)) {
250       if(fDebugMode)
251         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
252       return kFALSE;
253     }
254   }
255   if(fPointOnITSLayer5Flag) {
256     if(!track->HasPointOnITSLayer(4)) {
257       if(fDebugMode)
258         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
259       return kFALSE;
260     }
261   }
262   if(fPointOnITSLayer6Flag) {
263     if(!track->HasPointOnITSLayer(5)) {
264       if(fDebugMode)
265         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
266       return kFALSE;
267     }
268   }
269   if(fMinITSClustersFlag) {
270     if(nClustersITS < fMinITSClusters) {
271       if(fDebugMode)
272         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
273       return kFALSE;
274     }
275   }
276   if(fMaxChi2PerITSClusterFlag) {
277     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
278       if(fDebugMode)
279         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
280       return kFALSE; 
281     }
282   }
283   if(fMinTPCClustersFlag) {
284     if(nClustersTPC < fMinTPCClusters) {
285       if(fDebugMode)
286         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
287       return kFALSE;
288     }
289   }
290   if(fMaxChi2PerTPCClusterFlag) {
291     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
292       if(fDebugMode)
293         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
294       return kFALSE; 
295     }
296   }
297   if(fMaxCov11Flag) {
298     if(extCov[0] > fMaxCov11) {
299       if(fDebugMode)
300         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
301       return kFALSE;
302     }
303   }
304   if(fMaxCov22Flag) {
305     if(extCov[2] > fMaxCov22) {
306       if(fDebugMode)
307         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
308       return kFALSE;
309     }
310   }
311   if(fMaxCov33Flag) {
312     if(extCov[5] > fMaxCov33) {
313       if(fDebugMode)
314         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
315       return kFALSE;
316     }
317   }
318   if(fMaxCov44Flag) {
319     if(extCov[9] > fMaxCov44) {
320       if(fDebugMode)
321         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
322       return kFALSE;
323     }
324   }
325   if(fMaxCov55Flag) {
326     if(extCov[14] > fMaxCov55) {
327       if(fDebugMode)
328         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
329       return kFALSE;
330     }
331   }
332   if(fMaxSigmaToVertexFlag) {
333     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
334       if(fDebugMode)
335         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
336       return kFALSE;
337     }
338   }
339   if(fMaxSigmaToVertexTPCFlag) {
340     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
341       if(fDebugMode)
342         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
343       return kFALSE;
344     }
345   }
346   if(fMaxDCAXYFlag) { 
347     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
348       if(fDebugMode)
349         Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
350       return kFALSE;
351     }
352   }
353   if(fMaxDCAXYTPCFlag) { 
354     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
355       if(fDebugMode)
356         Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
357       return kFALSE;
358     }
359   }
360   if(fMaxDCAZFlag) { 
361     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
362       if(fDebugMode)
363         Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
364       return kFALSE;
365     }
366   }
367   if(fMaxDCAZTPCFlag) { 
368     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
369       if(fDebugMode)
370         Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
371       return kFALSE;
372     }
373   }
374   if(fMaxDCA3DFlag) { 
375     if(TMath::Abs(dca3D) > fMaxDCA3D) {
376       if(fDebugMode)
377         Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
378       return kFALSE;
379     }
380   }
381   if(fMaxDCA3DTPCFlag) { 
382     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
383       if(fDebugMode)
384         Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
385       return kFALSE;
386     }
387   }
388   if(fMaxConstrainChi2Flag) {
389     if(track->GetConstrainedChi2() > 0) 
390       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
391       if(fDebugMode)
392         Printf("IsAccepted: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
393       return kFALSE;
394       }
395   }
396   if(fMinTPCdEdxPointsFlag) {
397     if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
398       if(fDebugMode)
399         Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
400       return kFALSE;
401     }
402   }
403   if(fITSRefitFlag) {
404     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
405       if(fDebugMode)
406         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
407       return kFALSE;
408     }
409   }
410   if(fTPCRefitFlag) {
411     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
412       if(fDebugMode)
413         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
414       return kFALSE;
415     }
416   }
417   if(fESDpidFlag) {
418     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
419       if(fDebugMode)
420         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
421       return kFALSE;
422     }
423   }
424   if(fTPCpidFlag) {
425     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
426       if(fDebugMode)
427         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
428       return kFALSE;
429     }
430   }
431   if(fTOFpidFlag) {
432     if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
433       if(fDebugMode)
434         Printf("IsAccepted: Track rejected because it has no TOF pid flag");
435       return kFALSE;
436     }
437   }
438
439   return kTRUE;
440 }
441
442 //____________________________________________________________________//
443 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
444   // Calculates the number of sigma to the vertex.
445   
446   Float_t b[2];
447   Float_t bRes[2];
448   Float_t bCov[3];
449   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
450     esdTrack->GetImpactParametersTPC(b,bCov);
451   else
452     esdTrack->GetImpactParameters(b,bCov);
453   
454   if (bCov[0]<=0 || bCov[2]<=0) {
455     //AliDebug(1, "Estimated b resolution lower or equal zero!");
456     bCov[0]=0; bCov[2]=0;
457   }
458   bRes[0] = TMath::Sqrt(bCov[0]);
459   bRes[1] = TMath::Sqrt(bCov[2]);
460   
461   if (bRes[0] == 0 || bRes[1] ==0) return -1;
462   
463   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
464   
465   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
466   
467   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
468   
469   return d;
470 }
471
472 //____________________________________________________________________//
473 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx, 
474                                          Double_t gPy, 
475                                          Double_t gPz) const {
476   //returns the rapidity of the proton - to be removed
477   Double_t fMass = 9.38270000000000048e-01;
478   
479   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) + 
480                            TMath::Power(gPy,2) + 
481                            TMath::Power(gPz,2));
482   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
483   Double_t y = -999;
484   if(energy != gPz) 
485     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
486
487   return y;
488 }
489
490 //________________________________________________________________________
491 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd, 
492                                                      AnalysisMode mode,
493                                                      Double_t gVxMax,
494                                                      Double_t gVyMax,
495                                                      Double_t gVzMax) {
496   // Get the vertex from the ESD and returns it if the vertex is valid
497   // Second argument decides which vertex is used (this selects
498   // also the quality criteria that are applied)
499   const AliESDVertex* vertex = 0;
500   if(mode == kHybrid)
501     vertex = esd->GetPrimaryVertexSPD();
502   else if(mode == kTPC){
503     Double_t kBz = esd->GetMagneticField();
504     AliVertexerTracks vertexer(kBz);
505     vertexer.SetTPCMode();
506     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
507     esd->SetPrimaryVertexTPC(vTPC);
508     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
509       AliESDtrack *t = esd->GetTrack(i);
510       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
511     }
512     delete vTPC;
513     vertex = esd->GetPrimaryVertexTPC();
514   }
515   else if(mode == kGlobal)
516     vertex = esd->GetPrimaryVertex();
517   else
518     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
519   
520   if(!vertex) {
521     if(fDebugMode)
522       Printf("GetVertex: Event rejected because there is no valid vertex object");
523     return 0;
524   }
525   
526   // check Ncontributors
527   if(vertex->GetNContributors() <= 0) {
528     if(fDebugMode)
529       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
530     return 0;
531   }
532   
533   // check resolution
534   Double_t zRes = vertex->GetZRes();
535   if(zRes == 0) {
536     if(fDebugMode)
537       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
538     return 0;
539   }
540   
541   //check position
542   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
543     if(fDebugMode)
544       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
545     return 0;
546   }
547   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
548     if(fDebugMode)
549       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
550     return 0;
551   }
552   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
553     if(fDebugMode)
554       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
555     return 0;
556   }
557   
558   return vertex;
559 }
560
561 //________________________________________________________________________
562 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd, 
563                                                TriggerMode trigger) {
564   // check if the event was triggered
565   ULong64_t triggerMask = esd->GetTriggerMask();
566   
567   // definitions from p-p.cfg
568   ULong64_t spdFO = (1 << 14);
569   ULong64_t v0left = (1 << 11);
570   ULong64_t v0right = (1 << 12);
571
572   switch (trigger) {
573   case kMB1: {
574     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
575       return kTRUE;
576     break;
577   }
578   case kMB2: {
579     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
580       return kTRUE;
581     break;
582   }
583   case kSPDFASTOR: {
584     if (triggerMask & spdFO)
585       return kTRUE;
586     break;
587   }
588   }//switch
589
590   return kFALSE;
591 }
592
593 //________________________________________________________________________
594 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
595   // return the list of cuts and their cut values for reference
596   TLatex l;
597   l.SetTextAlign(12);
598   l.SetTextSize(0.04);
599
600   TString listOfCuts;
601
602   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
603   c->SetFillColor(10); c->SetHighLightColor(41);
604   c->Divide(3,2);
605
606   c->cd(1)->SetFillColor(10);
607   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
608
609   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
610   l.DrawLatex(0.1,0.82,listOfCuts.Data());
611   listOfCuts = "Analysis mode: ";
612   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone"; 
613   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC"; 
614   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking"; 
615   l.DrawLatex(0.1,0.74,listOfCuts.Data());
616   listOfCuts = "Trigger mode: "; 
617   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1"; 
618   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2"; 
619   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)"; 
620   l.DrawLatex(0.1,0.66,listOfCuts.Data());
621   listOfCuts = "PID mode: "; 
622   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
623   if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})"; 
624   if(fProtonPIDMode == kSigma1) {
625     listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
626     listOfCuts += " #sigma";
627   }
628   if(fProtonPIDMode == kSigma2) {
629     listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
630     listOfCuts += " #sigma";
631   }
632   l.DrawLatex(0.1,0.58,listOfCuts.Data());
633   listOfCuts = "Accepted vertex diamond: "; 
634   l.DrawLatex(0.1,0.5,listOfCuts.Data());
635   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
636   l.DrawLatex(0.6,0.5,listOfCuts.Data());
637   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
638   l.DrawLatex(0.6,0.4,listOfCuts.Data());
639   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
640   l.DrawLatex(0.6,0.3,listOfCuts.Data());
641   listOfCuts = "Phase space: "; 
642   l.DrawLatex(0.1,0.2,listOfCuts.Data());
643   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";  
644   else listOfCuts = "|y| < ";
645   listOfCuts += TMath::Abs(fMinX); 
646   l.DrawLatex(0.35,0.2,listOfCuts.Data());
647   listOfCuts = "N_{bins} = ";
648   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
649   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
650   l.DrawLatex(0.35,0.15,listOfCuts.Data());
651   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
652   listOfCuts += fMaxY; listOfCuts += "GeV/c"; 
653   l.DrawLatex(0.35,0.1,listOfCuts.Data());
654   listOfCuts = "N_{bins} = ";
655   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
656   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
657   l.DrawLatex(0.35,0.05,listOfCuts.Data());
658
659   c->cd(2)->SetFillColor(2);
660   l.DrawLatex(0.3,0.9,"ITS related cuts");
661   listOfCuts = "Request a cluster on SPD1: "; 
662   listOfCuts += fPointOnITSLayer1Flag;
663   l.DrawLatex(0.1,0.8,listOfCuts.Data());
664   listOfCuts = "Request a cluster on SPD2: "; 
665   listOfCuts += fPointOnITSLayer2Flag;
666   l.DrawLatex(0.1,0.7,listOfCuts.Data());
667   listOfCuts = "Request a cluster on SDD1: "; 
668   listOfCuts += fPointOnITSLayer3Flag;
669   l.DrawLatex(0.1,0.6,listOfCuts.Data());
670   listOfCuts = "Request a cluster on SDD2: "; 
671   listOfCuts += fPointOnITSLayer4Flag;
672   l.DrawLatex(0.1,0.5,listOfCuts.Data());
673   listOfCuts = "Request a cluster on SSD1: "; 
674   listOfCuts += fPointOnITSLayer5Flag;
675   l.DrawLatex(0.1,0.4,listOfCuts.Data());
676   listOfCuts = "Request a cluster on SSD2: "; 
677   listOfCuts += fPointOnITSLayer6Flag; 
678   l.DrawLatex(0.1,0.3,listOfCuts.Data());  
679   listOfCuts = "Minimum number of ITS clusters: ";
680   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
681   else listOfCuts += "Not used";
682   l.DrawLatex(0.1,0.2,listOfCuts.Data());
683   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
684   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster; 
685   else listOfCuts += "Not used";
686   l.DrawLatex(0.1,0.1,listOfCuts.Data());
687
688   c->cd(3)->SetFillColor(3);
689   l.DrawLatex(0.3,0.9,"TPC related cuts");
690   listOfCuts = "Minimum number of TPC clusters: ";
691   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
692   else listOfCuts += "Not used";
693   l.DrawLatex(0.1,0.7,listOfCuts.Data());
694   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
695   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
696   else listOfCuts += "Not used";
697   l.DrawLatex(0.1,0.5,listOfCuts.Data());
698   listOfCuts = "Minimum number of TPC points for the dE/dx: ";
699   if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
700   else listOfCuts += "Not used";
701   l.DrawLatex(0.1,0.3,listOfCuts.Data());
702
703   c->cd(4)->SetFillColor(4);
704   l.DrawLatex(0.3,0.9,"Tracking related cuts");
705   listOfCuts = "Maximum cov11: ";
706   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
707   else listOfCuts += "Not used";
708   l.DrawLatex(0.1,0.75,listOfCuts.Data());
709   listOfCuts = "Maximum cov22: ";
710   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
711   else listOfCuts += "Not used";
712   l.DrawLatex(0.1,0.6,listOfCuts.Data());
713   listOfCuts = "Maximum cov33: ";
714   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
715   else listOfCuts += "Not used";
716   l.DrawLatex(0.1,0.45,listOfCuts.Data());
717   listOfCuts = "Maximum cov44: ";
718   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
719   else listOfCuts += "Not used";
720   l.DrawLatex(0.1,0.3,listOfCuts.Data());
721   listOfCuts = "Maximum cov55: ";
722   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
723   else listOfCuts += "Not used";
724   l.DrawLatex(0.1,0.15,listOfCuts.Data());
725
726   c->cd(5)->SetFillColor(5);
727   l.DrawLatex(0.3,0.9,"DCA related cuts");
728   listOfCuts = "Maximum sigma to vertex: ";
729   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
730   else listOfCuts += "Not used";
731   l.DrawLatex(0.1,0.8,listOfCuts.Data());
732   listOfCuts = "Maximum sigma to vertex (TPC): ";
733   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
734   else listOfCuts += "Not used";
735   l.DrawLatex(0.1,0.72,listOfCuts.Data());
736   listOfCuts = "Maximum DCA in xy: ";
737   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
738   else listOfCuts += "Not used";
739   l.DrawLatex(0.1,0.64,listOfCuts.Data());
740   listOfCuts = "Maximum DCA in z: ";
741   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
742   else listOfCuts += "Not used";
743   l.DrawLatex(0.1,0.56,listOfCuts.Data());
744   listOfCuts = "Maximum DCA in 3D: ";
745   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
746   else listOfCuts += "Not used";
747   l.DrawLatex(0.1,0.48,listOfCuts.Data());
748   listOfCuts = "Maximum DCA in xy (TPC): ";
749   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
750   else listOfCuts += "Not used";
751   l.DrawLatex(0.1,0.4,listOfCuts.Data());
752   listOfCuts = "Maximum DCA in z (TPC): ";
753   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
754   else listOfCuts += "Not used";
755   l.DrawLatex(0.1,0.32,listOfCuts.Data());
756   listOfCuts = "Maximum DCA in 3D (TPC): ";
757   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
758   else listOfCuts += "Not used";
759   l.DrawLatex(0.1,0.24,listOfCuts.Data());
760   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
761   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
762   else listOfCuts += "Not used";
763   l.DrawLatex(0.1,0.16,listOfCuts.Data());
764
765   c->cd(6)->SetFillColor(6);
766   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
767   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
768   l.DrawLatex(0.1,0.7,listOfCuts.Data());
769   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag; 
770   l.DrawLatex(0.1,0.5,listOfCuts.Data());
771   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
772   l.DrawLatex(0.1,0.3,listOfCuts.Data());
773   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
774   l.DrawLatex(0.1,0.1,listOfCuts.Data());
775
776   return c;
777 }
778
779 //________________________________________________________________________
780 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
781   //Function that checks if a track is a proton
782   Double_t probability[5];
783   Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
784   Long64_t fParticleType = 0;
785  
786   //Bayesian approach for the PID
787   if(fProtonPIDMode == kBayesian) {
788     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
789       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
790       if(tpcTrack) {
791         gPt = tpcTrack->Pt();
792         gP = tpcTrack->P();
793         track->GetTPCpid(probability);
794       }
795     }//TPC standalone or Hybrid TPC
796     else if(fProtonAnalysisMode == kGlobal) {
797       gPt = track->Pt();
798       gP = track->P();
799       track->GetESDpid(probability);
800     }//Global tracking    
801     
802     Double_t rcc = 0.0;
803     for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
804       rcc += probability[i]*GetParticleFraction(i,gP);
805     if(rcc != 0.0) {
806       Double_t w[5];
807       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
808         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
809       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
810     }
811     if(fParticleType == 4)
812       return kTRUE;
813   }
814   //Ratio of the measured over the theoretical dE/dx a la STAR
815   else if(fProtonPIDMode == kRatio) {
816     Printf("The kRatio mode is not implemented yet!!!");
817     return kFALSE;
818   }
819   //Definition of an N-sigma area around the dE/dx vs P band
820   else if(fProtonPIDMode == kSigma1) {
821     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
822     if(tpcTrack) {
823       gPt = tpcTrack->Pt();
824       gP = tpcTrack->P();
825       gEta = tpcTrack->Eta();
826     }
827     //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
828     Int_t nbinP = Int_t(gP/0.05) - 6;
829     Double_t tpcSignal = track->GetTPCsignal();
830     Double_t dEdxMean = fdEdxMean[nbinP];
831     Double_t dEdxSigma = fdEdxSigma[nbinP];
832     if((tpcSignal <= dEdxMean + fNSigma*dEdxSigma)&&
833        (tpcSignal >= dEdxMean - fNSigma*dEdxSigma))
834       return kTRUE;
835   }//kSigma1 PID method
836   //Another definition of an N-sigma area around the dE/dx vs P band
837   else if(fProtonPIDMode == kSigma2) {
838     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
839     if(tpcTrack) {
840       gPt = tpcTrack->Pt();
841       gP = tpcTrack->P();
842       gEta = tpcTrack->Eta();
843     }
844     //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
845     Int_t nbinP = Int_t(gP/0.05) - 6;
846     Double_t tpcSignal = track->GetTPCsignal();
847     Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
848     Double_t dEdxSigma = fdEdxSigma[nbinP];
849     Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
850     if(nsigma <= fNSigma) 
851       return kTRUE;
852   }
853
854   return kFALSE;
855 }
856
857 //________________________________________________________________________
858 void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
859   // This function is used in case the kSigma1 or kSigma2 PID mode is selected
860   // It takes as an argument the name of the ascii file (for the time being) 
861   // that is generated as a prior process.
862   // This ascii file has three columns: The min. P value (bins of 50MeV/c) 
863   // the mean and the sigma of the dE/dx distributions for protons coming 
864   // from a gaussian fit.
865   ifstream in;
866   in.open(filename);
867
868   Double_t gPtMin = 0.0;
869   Int_t iCounter = 0;
870   while(in.good()) {
871     in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
872     if(fDebugMode)
873       Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
874     iCounter += 1;
875   }
876 }
877
878 //________________________________________________________________________
879 Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
880   // This is the Bethe-Bloch function normalised to 1 at the minimum
881   // We renormalize it based on the MC information
882   // WARNING: To be revised soon!!!
883   // This is just a temporary fix
884   Double_t normalization = 49.2;
885   Double_t bg2=bg*bg;
886   Double_t beta2 = bg2/(1.+ bg2);
887   Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
888   //
889   const Float_t kmeanCorrection =0.1;
890   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
891   bb *= meanCorrection;
892
893   return normalization*bb; 
894 }
895
896