]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG2/SPECTRA/AliProtonAnalysisBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Bug fix in ITS standalone tracking when reusing all clusters
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysisBase.cxx
... / ...
CommitLineData
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>
33class AliLog;
34class AliESDVertex;
35
36#include "AliProtonAnalysisBase.h"
37
38ClassImp(AliProtonAnalysisBase)
39
40//____________________________________________________________________//
41AliProtonAnalysisBase::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//____________________________________________________________________//
122AliProtonAnalysisBase::~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//____________________________________________________________________//
133Double_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//____________________________________________________________________//
149Bool_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//____________________________________________________________________//
206Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
207 const AliESDVertex *vertex,
208 AliESDtrack* track) {
209 // Checks if the track is excluded from the cuts
210 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
211 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
212 Double_t dca3D = 0.0;
213
214 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
215 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
216 if(!tpcTrack) {
217 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
218 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
219 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
220 }
221 else {
222 gPt = tpcTrack->Pt();
223 gPx = tpcTrack->Px();
224 gPy = tpcTrack->Py();
225 gPz = tpcTrack->Pz();
226 tpcTrack->PropagateToDCA(vertex,
227 esd->GetMagneticField(),
228 100.,dca,cov);
229 }
230 }//standalone TPC or hybrid TPC approaches
231 else {
232 gPt = track->Pt();
233 gPx = track->Px();
234 gPy = track->Py();
235 gPz = track->Pz();
236 track->PropagateToDCA(vertex,
237 esd->GetMagneticField(),
238 100.,dca,cov);
239 }
240 dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
241 TMath::Power(dca[1],2));
242
243 Int_t fIdxInt[200];
244 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
245 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
246
247 Float_t chi2PerClusterITS = -1;
248 if (nClustersITS!=0)
249 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
250 Float_t chi2PerClusterTPC = -1;
251 if (nClustersTPC!=0)
252 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
253
254 Double_t extCov[15];
255 track->GetExternalCovariance(extCov);
256
257 if(fPointOnITSLayer1Flag) {
258 if(!track->HasPointOnITSLayer(0)) {
259 if(fDebugMode)
260 Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
261 return kFALSE;
262 }
263 }
264 if(fPointOnITSLayer2Flag) {
265 if(!track->HasPointOnITSLayer(1)) {
266 if(fDebugMode)
267 Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
268 return kFALSE;
269 }
270 }
271 if(fPointOnITSLayer3Flag) {
272 if(!track->HasPointOnITSLayer(2)) {
273 if(fDebugMode)
274 Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
275 return kFALSE;
276 }
277 }
278 if(fPointOnITSLayer4Flag) {
279 if(!track->HasPointOnITSLayer(3)) {
280 if(fDebugMode)
281 Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
282 return kFALSE;
283 }
284 }
285 if(fPointOnITSLayer5Flag) {
286 if(!track->HasPointOnITSLayer(4)) {
287 if(fDebugMode)
288 Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
289 return kFALSE;
290 }
291 }
292 if(fPointOnITSLayer6Flag) {
293 if(!track->HasPointOnITSLayer(5)) {
294 if(fDebugMode)
295 Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
296 return kFALSE;
297 }
298 }
299 if(fMinITSClustersFlag) {
300 if(nClustersITS < fMinITSClusters) {
301 if(fDebugMode)
302 Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
303 return kFALSE;
304 }
305 }
306 if(fMaxChi2PerITSClusterFlag) {
307 if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
308 if(fDebugMode)
309 Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
310 return kFALSE;
311 }
312 }
313 if(fMinTPCClustersFlag) {
314 if(nClustersTPC < fMinTPCClusters) {
315 if(fDebugMode)
316 Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
317 return kFALSE;
318 }
319 }
320 if(fMaxChi2PerTPCClusterFlag) {
321 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
322 if(fDebugMode)
323 Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
324 return kFALSE;
325 }
326 }
327 if(fMaxCov11Flag) {
328 if(extCov[0] > fMaxCov11) {
329 if(fDebugMode)
330 Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
331 return kFALSE;
332 }
333 }
334 if(fMaxCov22Flag) {
335 if(extCov[2] > fMaxCov22) {
336 if(fDebugMode)
337 Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
338 return kFALSE;
339 }
340 }
341 if(fMaxCov33Flag) {
342 if(extCov[5] > fMaxCov33) {
343 if(fDebugMode)
344 Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
345 return kFALSE;
346 }
347 }
348 if(fMaxCov44Flag) {
349 if(extCov[9] > fMaxCov44) {
350 if(fDebugMode)
351 Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
352 return kFALSE;
353 }
354 }
355 if(fMaxCov55Flag) {
356 if(extCov[14] > fMaxCov55) {
357 if(fDebugMode)
358 Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
359 return kFALSE;
360 }
361 }
362 if(fMaxSigmaToVertexFlag) {
363 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
364 if(fDebugMode)
365 Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
366 return kFALSE;
367 }
368 }
369 if(fMaxSigmaToVertexTPCFlag) {
370 if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
371 if(fDebugMode)
372 Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
373 return kFALSE;
374 }
375 }
376 if(fMaxDCAXYFlag) {
377 if(TMath::Abs(dca[0]) > fMaxDCAXY) {
378 if(fDebugMode)
379 Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
380 return kFALSE;
381 }
382 }
383 if(fMaxDCAXYTPCFlag) {
384 if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
385 if(fDebugMode)
386 Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
387 return kFALSE;
388 }
389 }
390 if(fMaxDCAZFlag) {
391 if(TMath::Abs(dca[1]) > fMaxDCAZ) {
392 if(fDebugMode)
393 Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
394 return kFALSE;
395 }
396 }
397 if(fMaxDCAZTPCFlag) {
398 if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
399 if(fDebugMode)
400 Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
401 return kFALSE;
402 }
403 }
404 if(fMaxDCA3DFlag) {
405 if(TMath::Abs(dca3D) > fMaxDCA3D) {
406 if(fDebugMode)
407 Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
408 return kFALSE;
409 }
410 }
411 if(fMaxDCA3DTPCFlag) {
412 if(TMath::Abs(dca3D) > fMaxDCA3DTPC) {
413 if(fDebugMode)
414 Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
415 return kFALSE;
416 }
417 }
418 if(fMaxConstrainChi2Flag) {
419 if(track->GetConstrainedChi2() > 0)
420 if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) {
421 if(fDebugMode)
422 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);
423 return kFALSE;
424 }
425 }
426 if(fMinTPCdEdxPointsFlag) {
427 if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
428 if(fDebugMode)
429 Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
430 return kFALSE;
431 }
432 }
433 if(fITSRefitFlag) {
434 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
435 if(fDebugMode)
436 Printf("IsAccepted: Track rejected because it has no ITS refit flag");
437 return kFALSE;
438 }
439 }
440 if(fTPCRefitFlag) {
441 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
442 if(fDebugMode)
443 Printf("IsAccepted: Track rejected because it has no TPC refit flag");
444 return kFALSE;
445 }
446 }
447 if(fESDpidFlag) {
448 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
449 if(fDebugMode)
450 Printf("IsAccepted: Track rejected because it has no ESD pid flag");
451 return kFALSE;
452 }
453 }
454 if(fTPCpidFlag) {
455 if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
456 if(fDebugMode)
457 Printf("IsAccepted: Track rejected because it has no TPC pid flag");
458 return kFALSE;
459 }
460 }
461 if(fTOFpidFlag) {
462 if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
463 if(fDebugMode)
464 Printf("IsAccepted: Track rejected because it has no TOF pid flag");
465 return kFALSE;
466 }
467 }
468
469 return kTRUE;
470}
471
472//____________________________________________________________________//
473Bool_t AliProtonAnalysisBase::IsPrimary(AliESDEvent *esd,
474 const AliESDVertex *vertex,
475 AliESDtrack* track) {
476 // Checks if the track is a primary-like candidate
477 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
478 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
479 Double_t dca3D = 0.0;
480
481 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
482 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
483 if(!tpcTrack) {
484 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
485 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
486 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
487 }
488 else {
489 gPt = tpcTrack->Pt();
490 gPx = tpcTrack->Px();
491 gPy = tpcTrack->Py();
492 gPz = tpcTrack->Pz();
493 tpcTrack->PropagateToDCA(vertex,
494 esd->GetMagneticField(),
495 100.,dca,cov);
496 }
497 }//standalone TPC or hybrid TPC approaches
498 else {
499 gPt = track->Pt();
500 gPx = track->Px();
501 gPy = track->Py();
502 gPz = track->Pz();
503 track->PropagateToDCA(vertex,
504 esd->GetMagneticField(),
505 100.,dca,cov);
506 }
507 dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
508 TMath::Power(dca[1],2));
509
510 Int_t fIdxInt[200];
511 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
512 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
513
514 Float_t chi2PerClusterITS = -1;
515 if (nClustersITS!=0)
516 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
517 Float_t chi2PerClusterTPC = -1;
518 if (nClustersTPC!=0)
519 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
520
521 Double_t extCov[15];
522 track->GetExternalCovariance(extCov);
523
524 if(fPointOnITSLayer1Flag) {
525 if(!track->HasPointOnITSLayer(0)) {
526 if(fDebugMode)
527 Printf("IsPrimary: Track rejected because it doesn't have a point on the 1st ITS layer");
528 return kFALSE;
529 }
530 }
531 if(fPointOnITSLayer2Flag) {
532 if(!track->HasPointOnITSLayer(1)) {
533 if(fDebugMode)
534 Printf("IsPrimary: Track rejected because it doesn't have a point on the 2nd ITS layer");
535 return kFALSE;
536 }
537 }
538 if(fPointOnITSLayer3Flag) {
539 if(!track->HasPointOnITSLayer(2)) {
540 if(fDebugMode)
541 Printf("IsPrimary: Track rejected because it doesn't have a point on the 3rd ITS layer");
542 return kFALSE;
543 }
544 }
545 if(fPointOnITSLayer4Flag) {
546 if(!track->HasPointOnITSLayer(3)) {
547 if(fDebugMode)
548 Printf("IsPrimary: Track rejected because it doesn't have a point on the 4th ITS layer");
549 return kFALSE;
550 }
551 }
552 if(fPointOnITSLayer5Flag) {
553 if(!track->HasPointOnITSLayer(4)) {
554 if(fDebugMode)
555 Printf("IsPrimary: Track rejected because it doesn't have a point on the 5th ITS layer");
556 return kFALSE;
557 }
558 }
559 if(fPointOnITSLayer6Flag) {
560 if(!track->HasPointOnITSLayer(5)) {
561 if(fDebugMode)
562 Printf("IsPrimary: Track rejected because it doesn't have a point on the 6th ITS layer");
563 return kFALSE;
564 }
565 }
566 if(fMinITSClustersFlag) {
567 if(nClustersITS < fMinITSClusters) {
568 if(fDebugMode)
569 Printf("IsPrimary: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
570 return kFALSE;
571 }
572 }
573 if(fMaxChi2PerITSClusterFlag) {
574 if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
575 if(fDebugMode)
576 Printf("IsPrimary: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
577 return kFALSE;
578 }
579 }
580 if(fMinTPCClustersFlag) {
581 if(nClustersTPC < fMinTPCClusters) {
582 if(fDebugMode)
583 Printf("IsPrimary: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
584 return kFALSE;
585 }
586 }
587 if(fMaxChi2PerTPCClusterFlag) {
588 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
589 if(fDebugMode)
590 Printf("IsPrimary: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
591 return kFALSE;
592 }
593 }
594 if(fMaxCov11Flag) {
595 if(extCov[0] > fMaxCov11) {
596 if(fDebugMode)
597 Printf("IsPrimary: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
598 return kFALSE;
599 }
600 }
601 if(fMaxCov22Flag) {
602 if(extCov[2] > fMaxCov22) {
603 if(fDebugMode)
604 Printf("IsPrimary: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
605 return kFALSE;
606 }
607 }
608 if(fMaxCov33Flag) {
609 if(extCov[5] > fMaxCov33) {
610 if(fDebugMode)
611 Printf("IsPrimary: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
612 return kFALSE;
613 }
614 }
615 if(fMaxCov44Flag) {
616 if(extCov[9] > fMaxCov44) {
617 if(fDebugMode)
618 Printf("IsPrimary: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
619 return kFALSE;
620 }
621 }
622 if(fMaxCov55Flag) {
623 if(extCov[14] > fMaxCov55) {
624 if(fDebugMode)
625 Printf("IsPrimary: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
626 return kFALSE;
627 }
628 }
629 if(fMaxSigmaToVertexFlag) {
630 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
631 if(fDebugMode)
632 Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
633 return kFALSE;
634 }
635 }
636 if(fMaxSigmaToVertexTPCFlag) {
637 if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
638 if(fDebugMode)
639 Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
640 return kFALSE;
641 }
642 }
643 if(fMaxDCAXYFlag) {
644 if(TMath::Abs(dca[0]) > fMaxDCAXY) {
645 if(fDebugMode)
646 Printf("IsPrimary: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
647 return kFALSE;
648 }
649 }
650 if(fMaxDCAXYTPCFlag) {
651 if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
652 if(fDebugMode)
653 Printf("IsPrimary: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
654 return kFALSE;
655 }
656 }
657 if(fMaxDCAZFlag) {
658 if(TMath::Abs(dca[1]) > fMaxDCAZ) {
659 if(fDebugMode)
660 Printf("IsPrimary: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
661 return kFALSE;
662 }
663 }
664 if(fMaxDCAZTPCFlag) {
665 if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
666 if(fDebugMode)
667 Printf("IsPrimary: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
668 return kFALSE;
669 }
670 }
671 if(fMaxDCA3DFlag) {
672 if(TMath::Abs(dca3D) > fMaxDCA3D) {
673 if(fDebugMode)
674 Printf("IsPrimary: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
675 return kFALSE;
676 }
677 }
678 if(fMaxDCA3DTPCFlag) {
679 if(TMath::Abs(dca3D) > fMaxDCA3DTPC) {
680 if(fDebugMode)
681 Printf("IsPrimary: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
682 return kFALSE;
683 }
684 }
685 if(fMaxConstrainChi2Flag) {
686 if(track->GetConstrainedChi2() > 0)
687 if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) {
688 if(fDebugMode)
689 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);
690 return kFALSE;
691 }
692 }
693 if(fMinTPCdEdxPointsFlag) {
694 if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
695 if(fDebugMode)
696 Printf("IsPrimary: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
697 return kFALSE;
698 }
699 }
700 if(fITSRefitFlag) {
701 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
702 if(fDebugMode)
703 Printf("IsPrimary: Track rejected because it has no ITS refit flag");
704 return kFALSE;
705 }
706 }
707 if(fTPCRefitFlag) {
708 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
709 if(fDebugMode)
710 Printf("IsPrimary: Track rejected because it has no TPC refit flag");
711 return kFALSE;
712 }
713 }
714 if(fESDpidFlag) {
715 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
716 if(fDebugMode)
717 Printf("IsPrimary: Track rejected because it has no ESD pid flag");
718 return kFALSE;
719 }
720 }
721 if(fTPCpidFlag) {
722 if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
723 if(fDebugMode)
724 Printf("IsPrimary: Track rejected because it has no TPC pid flag");
725 return kFALSE;
726 }
727 }
728 if(fTOFpidFlag) {
729 if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
730 if(fDebugMode)
731 Printf("IsPrimary: Track rejected because it has no TOF pid flag");
732 return kFALSE;
733 }
734 }
735
736 return kTRUE;
737}
738
739//____________________________________________________________________//
740Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
741 // Calculates the number of sigma to the vertex.
742 Float_t b[2];
743 Float_t bRes[2];
744 Float_t bCov[3];
745 if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
746 esdTrack->GetImpactParametersTPC(b,bCov);
747 else
748 esdTrack->GetImpactParameters(b,bCov);
749
750 if (bCov[0]<=0 || bCov[2]<=0) {
751 //AliDebug(1, "Estimated b resolution lower or equal zero!");
752 bCov[0]=0; bCov[2]=0;
753 }
754 bRes[0] = TMath::Sqrt(bCov[0]);
755 bRes[1] = TMath::Sqrt(bCov[2]);
756
757 if (bRes[0] == 0 || bRes[1] ==0) return -1;
758
759 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
760
761 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
762
763 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
764
765 return d;
766}
767
768//____________________________________________________________________//
769Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
770 Double_t gPy,
771 Double_t gPz) const {
772 //returns the rapidity of the proton - to be removed
773 Double_t fMass = 9.38270000000000048e-01;
774
775 Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
776 TMath::Power(gPy,2) +
777 TMath::Power(gPz,2));
778 Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
779 Double_t y = -999;
780 if(energy != gPz)
781 y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
782
783 return y;
784}
785
786//________________________________________________________________________
787const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
788 AnalysisMode mode,
789 Double_t gVxMax,
790 Double_t gVyMax,
791 Double_t gVzMax) {
792 // Get the vertex from the ESD and returns it if the vertex is valid
793 // Second argument decides which vertex is used (this selects
794 // also the quality criteria that are applied)
795 const AliESDVertex* vertex = 0;
796 if(mode == kHybrid)
797 vertex = esd->GetPrimaryVertexSPD();
798 else if(mode == kTPC){
799 Double_t kBz = esd->GetMagneticField();
800 AliVertexerTracks vertexer(kBz);
801 vertexer.SetTPCMode();
802 AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
803 esd->SetPrimaryVertexTPC(vTPC);
804 for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
805 AliESDtrack *t = esd->GetTrack(i);
806 t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
807 }
808 delete vTPC;
809 vertex = esd->GetPrimaryVertexTPC();
810 }
811 else if(mode == kGlobal)
812 vertex = esd->GetPrimaryVertex();
813 else
814 Printf("GetVertex: ERROR: Invalid second argument %d", mode);
815
816 if(!vertex) {
817 if(fDebugMode)
818 Printf("GetVertex: Event rejected because there is no valid vertex object");
819 return 0;
820 }
821
822 // check Ncontributors
823 if(vertex->GetNContributors() <= 0) {
824 if(fDebugMode)
825 Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
826 return 0;
827 }
828
829 // check resolution
830 Double_t zRes = vertex->GetZRes();
831 if(zRes == 0) {
832 if(fDebugMode)
833 Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
834 return 0;
835 }
836 ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
837 ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
838 ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
839
840 //check position
841 if(TMath::Abs(vertex->GetXv()) > gVxMax) {
842 if(fDebugMode)
843 Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
844 return 0;
845 }
846 if(TMath::Abs(vertex->GetYv()) > gVyMax) {
847 if(fDebugMode)
848 Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
849 return 0;
850 }
851 if(TMath::Abs(vertex->GetZv()) > gVzMax) {
852 if(fDebugMode)
853 Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
854 return 0;
855 }
856 ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
857 ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
858 ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
859 ((TH1F *)(fListVertexQA->At(6)))->Fill(vertex->GetNContributors());
860
861 //check number of contributors
862 if(fMinNumOfContributors > 0) {
863 if(fMinNumOfContributors > vertex->GetNContributors()) {
864 if(fDebugMode)
865 Printf("GetVertex: Event rejected because it has %d number of contributors (requested minimum: %d)",vertex->GetNContributors(),fMinNumOfContributors);
866
867 return 0;
868 }
869 }
870
871 return vertex;
872}
873
874//________________________________________________________________________
875Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
876 TriggerMode trigger) {
877 // check if the event was triggered
878 ULong64_t triggerMask = esd->GetTriggerMask();
879 TString firedTriggerClass = esd->GetFiredTriggerClasses();
880
881 if(fAnalysisMC) {
882 // definitions from p-p.cfg
883 ULong64_t spdFO = (1 << 14);
884 ULong64_t v0left = (1 << 11);
885 ULong64_t v0right = (1 << 12);
886
887 switch (trigger) {
888 case kMB1: {
889 if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
890 return kTRUE;
891 break;
892 }
893 case kMB2: {
894 if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
895 return kTRUE;
896 break;
897 }
898 case kSPDFASTOR: {
899 if (triggerMask & spdFO)
900 return kTRUE;
901 break;
902 }
903 }//switch
904 }
905 else {
906 if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
907 return kTRUE;
908 }
909
910 return kFALSE;
911}
912
913//________________________________________________________________________
914TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
915 // return the list of cuts and their cut values for reference
916 TLatex l;
917 l.SetTextAlign(12);
918 l.SetTextSize(0.04);
919
920 TString listOfCuts;
921
922 TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
923 c->SetFillColor(10); c->SetHighLightColor(41);
924 c->Divide(3,2);
925
926 c->cd(1)->SetFillColor(10);
927 l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
928
929 listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
930 l.DrawLatex(0.1,0.82,listOfCuts.Data());
931 listOfCuts = "Analysis mode: ";
932 if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
933 if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
934 if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
935 l.DrawLatex(0.1,0.74,listOfCuts.Data());
936 listOfCuts = "Trigger mode: ";
937 if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
938 if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
939 if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
940 l.DrawLatex(0.1,0.66,listOfCuts.Data());
941 listOfCuts = "PID mode: ";
942 if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
943 if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})";
944 if(fProtonPIDMode == kSigma1) {
945 listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
946 listOfCuts += " #sigma";
947 }
948 if(fProtonPIDMode == kSigma2) {
949 listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
950 listOfCuts += " #sigma";
951 }
952 l.DrawLatex(0.1,0.58,listOfCuts.Data());
953 listOfCuts = "Accepted vertex diamond: ";
954 l.DrawLatex(0.1,0.5,listOfCuts.Data());
955 listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
956 l.DrawLatex(0.6,0.5,listOfCuts.Data());
957 listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
958 l.DrawLatex(0.6,0.4,listOfCuts.Data());
959 listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
960 l.DrawLatex(0.6,0.3,listOfCuts.Data());
961 listOfCuts = "Phase space: ";
962 l.DrawLatex(0.1,0.2,listOfCuts.Data());
963 if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
964 else listOfCuts = "|y| < ";
965 listOfCuts += TMath::Abs(fMinX);
966 l.DrawLatex(0.35,0.2,listOfCuts.Data());
967 listOfCuts = "N_{bins} = ";
968 listOfCuts += fNBinsX; listOfCuts += " (binning: ";
969 listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
970 l.DrawLatex(0.35,0.15,listOfCuts.Data());
971 listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
972 listOfCuts += fMaxY; listOfCuts += "GeV/c";
973 l.DrawLatex(0.35,0.1,listOfCuts.Data());
974 listOfCuts = "N_{bins} = ";
975 listOfCuts += fNBinsY; listOfCuts += " (binning: ";
976 listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
977 l.DrawLatex(0.35,0.05,listOfCuts.Data());
978
979 c->cd(2)->SetFillColor(2);
980 l.DrawLatex(0.3,0.9,"ITS related cuts");
981 listOfCuts = "Request a cluster on SPD1: ";
982 listOfCuts += fPointOnITSLayer1Flag;
983 l.DrawLatex(0.1,0.8,listOfCuts.Data());
984 listOfCuts = "Request a cluster on SPD2: ";
985 listOfCuts += fPointOnITSLayer2Flag;
986 l.DrawLatex(0.1,0.7,listOfCuts.Data());
987 listOfCuts = "Request a cluster on SDD1: ";
988 listOfCuts += fPointOnITSLayer3Flag;
989 l.DrawLatex(0.1,0.6,listOfCuts.Data());
990 listOfCuts = "Request a cluster on SDD2: ";
991 listOfCuts += fPointOnITSLayer4Flag;
992 l.DrawLatex(0.1,0.5,listOfCuts.Data());
993 listOfCuts = "Request a cluster on SSD1: ";
994 listOfCuts += fPointOnITSLayer5Flag;
995 l.DrawLatex(0.1,0.4,listOfCuts.Data());
996 listOfCuts = "Request a cluster on SSD2: ";
997 listOfCuts += fPointOnITSLayer6Flag;
998 l.DrawLatex(0.1,0.3,listOfCuts.Data());
999 listOfCuts = "Minimum number of ITS clusters: ";
1000 if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
1001 else listOfCuts += "Not used";
1002 l.DrawLatex(0.1,0.2,listOfCuts.Data());
1003 listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
1004 if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
1005 else listOfCuts += "Not used";
1006 l.DrawLatex(0.1,0.1,listOfCuts.Data());
1007
1008 c->cd(3)->SetFillColor(3);
1009 l.DrawLatex(0.3,0.9,"TPC related cuts");
1010 listOfCuts = "Minimum number of TPC clusters: ";
1011 if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
1012 else listOfCuts += "Not used";
1013 l.DrawLatex(0.1,0.7,listOfCuts.Data());
1014 listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
1015 if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
1016 else listOfCuts += "Not used";
1017 l.DrawLatex(0.1,0.5,listOfCuts.Data());
1018 listOfCuts = "Minimum number of TPC points for the dE/dx: ";
1019 if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
1020 else listOfCuts += "Not used";
1021 l.DrawLatex(0.1,0.3,listOfCuts.Data());
1022
1023 c->cd(4)->SetFillColor(4);
1024 l.DrawLatex(0.3,0.9,"Tracking related cuts");
1025 listOfCuts = "Maximum cov11: ";
1026 if(fMaxCov11Flag) listOfCuts += fMaxCov11;
1027 else listOfCuts += "Not used";
1028 l.DrawLatex(0.1,0.75,listOfCuts.Data());
1029 listOfCuts = "Maximum cov22: ";
1030 if(fMaxCov22Flag) listOfCuts += fMaxCov22;
1031 else listOfCuts += "Not used";
1032 l.DrawLatex(0.1,0.6,listOfCuts.Data());
1033 listOfCuts = "Maximum cov33: ";
1034 if(fMaxCov33Flag) listOfCuts += fMaxCov33;
1035 else listOfCuts += "Not used";
1036 l.DrawLatex(0.1,0.45,listOfCuts.Data());
1037 listOfCuts = "Maximum cov44: ";
1038 if(fMaxCov44Flag) listOfCuts += fMaxCov44;
1039 else listOfCuts += "Not used";
1040 l.DrawLatex(0.1,0.3,listOfCuts.Data());
1041 listOfCuts = "Maximum cov55: ";
1042 if(fMaxCov55Flag) listOfCuts += fMaxCov55;
1043 else listOfCuts += "Not used";
1044 l.DrawLatex(0.1,0.15,listOfCuts.Data());
1045
1046 c->cd(5)->SetFillColor(5);
1047 l.DrawLatex(0.3,0.9,"DCA related cuts");
1048 listOfCuts = "Maximum sigma to vertex: ";
1049 if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
1050 else listOfCuts += "Not used";
1051 l.DrawLatex(0.1,0.8,listOfCuts.Data());
1052 listOfCuts = "Maximum sigma to vertex (TPC): ";
1053 if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
1054 else listOfCuts += "Not used";
1055 l.DrawLatex(0.1,0.72,listOfCuts.Data());
1056 listOfCuts = "Maximum DCA in xy: ";
1057 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
1058 else listOfCuts += "Not used";
1059 l.DrawLatex(0.1,0.64,listOfCuts.Data());
1060 listOfCuts = "Maximum DCA in z: ";
1061 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
1062 else listOfCuts += "Not used";
1063 l.DrawLatex(0.1,0.56,listOfCuts.Data());
1064 listOfCuts = "Maximum DCA in 3D: ";
1065 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
1066 else listOfCuts += "Not used";
1067 l.DrawLatex(0.1,0.48,listOfCuts.Data());
1068 listOfCuts = "Maximum DCA in xy (TPC): ";
1069 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
1070 else listOfCuts += "Not used";
1071 l.DrawLatex(0.1,0.4,listOfCuts.Data());
1072 listOfCuts = "Maximum DCA in z (TPC): ";
1073 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
1074 else listOfCuts += "Not used";
1075 l.DrawLatex(0.1,0.32,listOfCuts.Data());
1076 listOfCuts = "Maximum DCA in 3D (TPC): ";
1077 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
1078 else listOfCuts += "Not used";
1079 l.DrawLatex(0.1,0.24,listOfCuts.Data());
1080 listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
1081 if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
1082 else listOfCuts += "Not used";
1083 l.DrawLatex(0.1,0.16,listOfCuts.Data());
1084
1085 c->cd(6)->SetFillColor(6);
1086 l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
1087 listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
1088 l.DrawLatex(0.1,0.7,listOfCuts.Data());
1089 listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
1090 l.DrawLatex(0.1,0.5,listOfCuts.Data());
1091 listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
1092 l.DrawLatex(0.1,0.3,listOfCuts.Data());
1093 listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
1094 l.DrawLatex(0.1,0.1,listOfCuts.Data());
1095
1096 return c;
1097}
1098
1099//________________________________________________________________________
1100Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
1101 //Function that checks if a track is a proton
1102 Double_t probability[5];
1103 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
1104 Long64_t fParticleType = 0;
1105
1106 //Bayesian approach for the PID
1107 if(fProtonPIDMode == kBayesian) {
1108 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
1109 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
1110 if(tpcTrack) {
1111 gPt = tpcTrack->Pt();
1112 gP = tpcTrack->P();
1113 track->GetTPCpid(probability);
1114 }
1115 }//TPC standalone or Hybrid TPC
1116 else if(fProtonAnalysisMode == kGlobal) {
1117 gPt = track->Pt();
1118 gP = track->P();
1119 track->GetESDpid(probability);
1120 }//Global tracking
1121
1122 Double_t rcc = 0.0;
1123 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
1124 rcc += probability[i]*GetParticleFraction(i,gP);
1125 if(rcc != 0.0) {
1126 Double_t w[5];
1127 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
1128 w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
1129 fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
1130 }
1131 if(fParticleType == 4)
1132 return kTRUE;
1133 }
1134 //Ratio of the measured over the theoretical dE/dx a la STAR
1135 else if(fProtonPIDMode == kRatio) {
1136 Printf("The kRatio mode is not implemented yet!!!");
1137 return kFALSE;
1138 }
1139 //Definition of an N-sigma area around the dE/dx vs P band
1140 else if(fProtonPIDMode == kSigma1) {
1141 Double_t fAlephParameters[5];
1142 if(fAnalysisMC) {
1143 fAlephParameters[0] = 2.15898e+00/50.;
1144 fAlephParameters[1] = 1.75295e+01;
1145 fAlephParameters[2] = 3.40030e-09;
1146 fAlephParameters[3] = 1.96178e+00;
1147 fAlephParameters[4] = 3.91720e+00;
1148 }
1149 else {
1150 fAlephParameters[0] = 0.0283086;
1151 fAlephParameters[1] = 2.63394e+01;
1152 fAlephParameters[2] = 5.04114e-11;
1153 fAlephParameters[3] = 2.12543e+00;
1154 fAlephParameters[4] = 4.88663e+00;
1155 }
1156
1157 AliESDpid *fESDpid = new AliESDpid();
1158 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
1159
1160 Double_t nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton));
1161 if(nsigma <= fNSigma)
1162 return kTRUE;
1163 }//kSigma1 PID method
1164 //Another definition of an N-sigma area around the dE/dx vs P band
1165 else if(fProtonPIDMode == kSigma2) {
1166 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
1167 if(tpcTrack) {
1168 gPt = tpcTrack->Pt();
1169 gP = tpcTrack->P();
1170 gEta = tpcTrack->Eta();
1171 }
1172 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
1173 Int_t nbinP = Int_t(gP/0.05) - 6;
1174 Double_t tpcSignal = track->GetTPCsignal();
1175 Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
1176 Double_t dEdxSigma = fdEdxSigma[nbinP];
1177 Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
1178 if(nsigma <= fNSigma)
1179 return kTRUE;
1180 }
1181
1182 return kFALSE;
1183}
1184
1185//________________________________________________________________________
1186void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
1187 // This function is used in case the kSigma1 or kSigma2 PID mode is selected
1188 // It takes as an argument the name of the ascii file (for the time being)
1189 // that is generated as a prior process.
1190 // This ascii file has three columns: The min. P value (bins of 50MeV/c)
1191 // the mean and the sigma of the dE/dx distributions for protons coming
1192 // from a gaussian fit.
1193 ifstream in;
1194 in.open(filename);
1195
1196 Double_t gPtMin = 0.0;
1197 Int_t iCounter = 0;
1198 while(in.good()) {
1199 in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
1200 if(fDebugMode)
1201 Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
1202 iCounter += 1;
1203 }
1204}
1205
1206//________________________________________________________________________
1207Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
1208 // This is the Bethe-Bloch function normalised to 1 at the minimum
1209 // We renormalize it based on the MC information
1210 // WARNING: To be revised soon!!!
1211 // This is just a temporary fix
1212 Double_t normalization = 49.2;
1213 Double_t bg2=bg*bg;
1214 Double_t beta2 = bg2/(1.+ bg2);
1215 Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
1216 //
1217 const Float_t kmeanCorrection =0.1;
1218 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
1219 bb *= meanCorrection;
1220
1221 return normalization*bb;
1222}
1223
1224