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