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