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