]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/AliProtonAnalysisBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix in ITS standalone tracking when reusing all clusters
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysisBase.cxx
CommitLineData
0ab648ea 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>
0ab648ea 22#include <TCanvas.h>
23#include <TLatex.h>
73aba974 24#include <TF1.h>
735cc63d 25#include <TList.h>
26#include <TH1F.h>
0ab648ea 27
28#include <AliExternalTrackParam.h>
29#include <AliESDEvent.h>
0ab648ea 30#include <AliPID.h>
0ab648ea 31#include <AliVertexerTracks.h>
10d100d4 32#include <AliESDpid.h>
73aba974 33class AliLog;
34class AliESDVertex;
35
36#include "AliProtonAnalysisBase.h"
37
0ab648ea 38ClassImp(AliProtonAnalysisBase)
39
40//____________________________________________________________________//
41AliProtonAnalysisBase::AliProtonAnalysisBase() :
790140ac 42 TObject(), fProtonAnalysisLevel("ESD"), fAnalysisMC(kFALSE),
e56f08ed 43 fTriggerMode(kMB2), kUseOnlineTrigger(kFALSE), kUseOfflineTrigger(kFALSE),
44 fPhysicsSelection(0),
76161fe0 45 fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
0ab648ea 46 fAnalysisEtaMode(kFALSE),
e56f08ed 47 fVxMax(100.), fVyMax(100.), fVzMax(100.), fMinNumOfContributors(0),
0ab648ea 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),
87a55728 57 fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
0ab648ea 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),
f62e9410 68 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
0ab648ea 69 fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
70 fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
71 fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
87a55728 72 fMinTPCdEdxPointsFlag(kFALSE),
0ab648ea 73 fFunctionProbabilityFlag(kFALSE),
87a55728 74 fNSigma(0),
0ab648ea 75 fElectronFunction(0), fMuonFunction(0),
76 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
df201289 77 fDebugMode(kFALSE), fListVertexQA(new TList()) {
0ab648ea 78 //Default constructor
79 for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
d488527a 80 for(Int_t i = 0; i < 24; i++) {
87a55728 81 fdEdxMean[i] = 0.0;
82 fdEdxSigma[i] = 0.0;
83 }
735cc63d 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
e56f08ed 113 TH1F *gHistNumberOfContributors = new TH1F("gHistNumberOfContributors",
114 "Number of contributors;N_{contr.};Entries",
115 100,0.,100.);
116 fListVertexQA->Add(gHistNumberOfContributors);
117
735cc63d 118
0ab648ea 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;
735cc63d 129 if(fListVertexQA) delete fListVertexQA;
0ab648ea 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
df201289 151 Double_t gP = 0.0, gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
0ab648ea 152 Double_t eta = 0.0;
153
a84984fa 154 if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid)) {
0ab648ea 155 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
156 if(!tpcTrack) {
df201289 157 gP = 0.0; gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
0ab648ea 158 }
159 else {
df201289 160 gP = tpcTrack->P();
0ab648ea 161 gPt = tpcTrack->Pt();
162 gPx = tpcTrack->Px();
163 gPy = tpcTrack->Py();
164 gPz = tpcTrack->Pz();
165 eta = tpcTrack->Eta();
166 }
a84984fa 167 }//standalone TPC or Hybrid TPC approaches
0ab648ea 168 else {
df201289 169 gP = track->P();
0ab648ea 170 gPt = track->Pt();
171 gPx = track->Px();
172 gPy = track->Py();
173 gPz = track->Pz();
174 eta = track->Eta();
175 }
176
e56f08ed 177 if((gPt < fMinY) || (gPt > fMaxY)) {
0ab648ea 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 }
e56f08ed 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 }
0ab648ea 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
a84984fa 214 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
0ab648ea 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 }
a84984fa 230 }//standalone TPC or hybrid TPC approaches
231 else {
0ab648ea 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 }
87a55728 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 }
0ab648ea 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 }
f62e9410 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 }
0ab648ea 468
469 return kTRUE;
470}
471
e56f08ed 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
0ab648ea 739//____________________________________________________________________//
740Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
741 // Calculates the number of sigma to the vertex.
0ab648ea 742 Float_t b[2];
743 Float_t bRes[2];
744 Float_t bCov[3];
ec785069 745 if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
0ab648ea 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 }
735cc63d 836 ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
837 ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
838 ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
e56f08ed 839
0ab648ea 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 }
735cc63d 856 ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
857 ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
858 ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
e56f08ed 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 }
0ab648ea 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();
42270c4c 879 TString firedTriggerClass = esd->GetFiredTriggerClasses();
0ab648ea 880
42270c4c 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
0ab648ea 904 }
42270c4c 905 else {
906 if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
0ab648ea 907 return kTRUE;
0ab648ea 908 }
0ab648ea 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.})";
87a55728 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 }
0ab648ea 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());
57e749bb 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());
0ab648ea 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];
87a55728 1103 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
0ab648ea 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
87a55728 1140 else if(fProtonPIDMode == kSigma1) {
2c080079 1141 Double_t fAlephParameters[5];
0bb8f6b4 1142 if(fAnalysisMC) {
e56f08ed 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;
0bb8f6b4 1148 }
1149 else {
e56f08ed 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;
0bb8f6b4 1155 }
2c080079 1156
10d100d4 1157 AliESDpid *fESDpid = new AliESDpid();
e56f08ed 1158 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
2c080079 1159
10d100d4 1160 Double_t nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton));
2c080079 1161 if(nsigma <= fNSigma)
87a55728 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;
0ab648ea 1180 }
1181
1182 return kFALSE;
1183}
1184
87a55728 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//________________________________________________________________________
73aba974 1207Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
87a55728 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}
0ab648ea 1223
1224