]> 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
old code removed
[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),
64 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
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;
87a55728 76 for(Int_t i = 0; i < 5; i++) {
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
114 if(kTPC) {
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
167 if((kTPC)&&(!kHybrid)) {
168 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
169 if(!tpcTrack) {
170 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
171 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
172 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
173 }
174 else {
175 gPt = tpcTrack->Pt();
176 gPx = tpcTrack->Px();
177 gPy = tpcTrack->Py();
178 gPz = tpcTrack->Pz();
179 tpcTrack->PropagateToDCA(vertex,
180 esd->GetMagneticField(),
181 100.,dca,cov);
182 }
183 }
184 else if(kHybrid) {
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 }
431
432 return kTRUE;
433}
434
435//____________________________________________________________________//
436Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
437 // Calculates the number of sigma to the vertex.
438
439 Float_t b[2];
440 Float_t bRes[2];
441 Float_t bCov[3];
442 if((kTPC)&&(!kHybrid))
443 esdTrack->GetImpactParametersTPC(b,bCov);
444 else
445 esdTrack->GetImpactParameters(b,bCov);
446
447 if (bCov[0]<=0 || bCov[2]<=0) {
448 //AliDebug(1, "Estimated b resolution lower or equal zero!");
449 bCov[0]=0; bCov[2]=0;
450 }
451 bRes[0] = TMath::Sqrt(bCov[0]);
452 bRes[1] = TMath::Sqrt(bCov[2]);
453
454 if (bRes[0] == 0 || bRes[1] ==0) return -1;
455
456 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
457
458 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
459
460 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
461
462 return d;
463}
464
465//____________________________________________________________________//
466Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
467 Double_t gPy,
468 Double_t gPz) const {
469 //returns the rapidity of the proton - to be removed
470 Double_t fMass = 9.38270000000000048e-01;
471
472 Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
473 TMath::Power(gPy,2) +
474 TMath::Power(gPz,2));
475 Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
476 Double_t y = -999;
477 if(energy != gPz)
478 y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
479
480 return y;
481}
482
483//________________________________________________________________________
484const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
485 AnalysisMode mode,
486 Double_t gVxMax,
487 Double_t gVyMax,
488 Double_t gVzMax) {
489 // Get the vertex from the ESD and returns it if the vertex is valid
490 // Second argument decides which vertex is used (this selects
491 // also the quality criteria that are applied)
492 const AliESDVertex* vertex = 0;
493 if(mode == kHybrid)
494 vertex = esd->GetPrimaryVertexSPD();
495 else if(mode == kTPC){
496 Double_t kBz = esd->GetMagneticField();
497 AliVertexerTracks vertexer(kBz);
498 vertexer.SetTPCMode();
499 AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
500 esd->SetPrimaryVertexTPC(vTPC);
501 for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
502 AliESDtrack *t = esd->GetTrack(i);
503 t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
504 }
505 delete vTPC;
506 vertex = esd->GetPrimaryVertexTPC();
507 }
508 else if(mode == kGlobal)
509 vertex = esd->GetPrimaryVertex();
510 else
511 Printf("GetVertex: ERROR: Invalid second argument %d", mode);
512
513 if(!vertex) {
514 if(fDebugMode)
515 Printf("GetVertex: Event rejected because there is no valid vertex object");
516 return 0;
517 }
518
519 // check Ncontributors
520 if(vertex->GetNContributors() <= 0) {
521 if(fDebugMode)
522 Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
523 return 0;
524 }
525
526 // check resolution
527 Double_t zRes = vertex->GetZRes();
528 if(zRes == 0) {
529 if(fDebugMode)
530 Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
531 return 0;
532 }
533
534 //check position
535 if(TMath::Abs(vertex->GetXv()) > gVxMax) {
536 if(fDebugMode)
537 Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
538 return 0;
539 }
540 if(TMath::Abs(vertex->GetYv()) > gVyMax) {
541 if(fDebugMode)
542 Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
543 return 0;
544 }
545 if(TMath::Abs(vertex->GetZv()) > gVzMax) {
546 if(fDebugMode)
547 Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
548 return 0;
549 }
550
551 return vertex;
552}
553
554//________________________________________________________________________
555Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
556 TriggerMode trigger) {
557 // check if the event was triggered
558 ULong64_t triggerMask = esd->GetTriggerMask();
559
560 // definitions from p-p.cfg
561 ULong64_t spdFO = (1 << 14);
562 ULong64_t v0left = (1 << 11);
563 ULong64_t v0right = (1 << 12);
564
565 switch (trigger) {
566 case kMB1: {
567 if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
568 return kTRUE;
569 break;
570 }
571 case kMB2: {
572 if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
573 return kTRUE;
574 break;
575 }
576 case kSPDFASTOR: {
577 if (triggerMask & spdFO)
578 return kTRUE;
579 break;
580 }
581 }//switch
582
583 return kFALSE;
584}
585
586//________________________________________________________________________
587TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
588 // return the list of cuts and their cut values for reference
589 TLatex l;
590 l.SetTextAlign(12);
591 l.SetTextSize(0.04);
592
593 TString listOfCuts;
594
595 TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
596 c->SetFillColor(10); c->SetHighLightColor(41);
597 c->Divide(3,2);
598
599 c->cd(1)->SetFillColor(10);
600 l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
601
602 listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
603 l.DrawLatex(0.1,0.82,listOfCuts.Data());
604 listOfCuts = "Analysis mode: ";
605 if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
606 if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
607 if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
608 l.DrawLatex(0.1,0.74,listOfCuts.Data());
609 listOfCuts = "Trigger mode: ";
610 if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
611 if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
612 if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
613 l.DrawLatex(0.1,0.66,listOfCuts.Data());
614 listOfCuts = "PID mode: ";
615 if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
616 if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})";
87a55728 617 if(fProtonPIDMode == kSigma1) {
618 listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
619 listOfCuts += " #sigma";
620 }
621 if(fProtonPIDMode == kSigma2) {
622 listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
623 listOfCuts += " #sigma";
624 }
0ab648ea 625 l.DrawLatex(0.1,0.58,listOfCuts.Data());
626 listOfCuts = "Accepted vertex diamond: ";
627 l.DrawLatex(0.1,0.5,listOfCuts.Data());
628 listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
629 l.DrawLatex(0.6,0.5,listOfCuts.Data());
630 listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
631 l.DrawLatex(0.6,0.4,listOfCuts.Data());
632 listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
633 l.DrawLatex(0.6,0.3,listOfCuts.Data());
634 listOfCuts = "Phase space: ";
635 l.DrawLatex(0.1,0.2,listOfCuts.Data());
636 if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
637 else listOfCuts = "|y| < ";
638 listOfCuts += TMath::Abs(fMinX);
639 l.DrawLatex(0.35,0.2,listOfCuts.Data());
640 listOfCuts = "N_{bins} = ";
641 listOfCuts += fNBinsX; listOfCuts += " (binning: ";
642 listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
643 l.DrawLatex(0.35,0.15,listOfCuts.Data());
644 listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
645 listOfCuts += fMaxY; listOfCuts += "GeV/c";
646 l.DrawLatex(0.35,0.1,listOfCuts.Data());
647 listOfCuts = "N_{bins} = ";
648 listOfCuts += fNBinsY; listOfCuts += " (binning: ";
649 listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
650 l.DrawLatex(0.35,0.05,listOfCuts.Data());
651
652 c->cd(2)->SetFillColor(2);
653 l.DrawLatex(0.3,0.9,"ITS related cuts");
654 listOfCuts = "Request a cluster on SPD1: ";
655 listOfCuts += fPointOnITSLayer1Flag;
656 l.DrawLatex(0.1,0.8,listOfCuts.Data());
657 listOfCuts = "Request a cluster on SPD2: ";
658 listOfCuts += fPointOnITSLayer2Flag;
659 l.DrawLatex(0.1,0.7,listOfCuts.Data());
660 listOfCuts = "Request a cluster on SDD1: ";
661 listOfCuts += fPointOnITSLayer3Flag;
662 l.DrawLatex(0.1,0.6,listOfCuts.Data());
663 listOfCuts = "Request a cluster on SDD2: ";
664 listOfCuts += fPointOnITSLayer4Flag;
665 l.DrawLatex(0.1,0.5,listOfCuts.Data());
666 listOfCuts = "Request a cluster on SSD1: ";
667 listOfCuts += fPointOnITSLayer5Flag;
668 l.DrawLatex(0.1,0.4,listOfCuts.Data());
669 listOfCuts = "Request a cluster on SSD2: ";
670 listOfCuts += fPointOnITSLayer6Flag;
671 l.DrawLatex(0.1,0.3,listOfCuts.Data());
672 listOfCuts = "Minimum number of ITS clusters: ";
673 if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
674 else listOfCuts += "Not used";
675 l.DrawLatex(0.1,0.2,listOfCuts.Data());
676 listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
677 if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
678 else listOfCuts += "Not used";
679 l.DrawLatex(0.1,0.1,listOfCuts.Data());
680
681 c->cd(3)->SetFillColor(3);
682 l.DrawLatex(0.3,0.9,"TPC related cuts");
683 listOfCuts = "Minimum number of TPC clusters: ";
684 if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
685 else listOfCuts += "Not used";
686 l.DrawLatex(0.1,0.7,listOfCuts.Data());
687 listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
688 if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
689 else listOfCuts += "Not used";
690 l.DrawLatex(0.1,0.5,listOfCuts.Data());
691
692 c->cd(4)->SetFillColor(4);
693 l.DrawLatex(0.3,0.9,"Tracking related cuts");
694 listOfCuts = "Maximum cov11: ";
695 if(fMaxCov11Flag) listOfCuts += fMaxCov11;
696 else listOfCuts += "Not used";
697 l.DrawLatex(0.1,0.75,listOfCuts.Data());
698 listOfCuts = "Maximum cov22: ";
699 if(fMaxCov22Flag) listOfCuts += fMaxCov22;
700 else listOfCuts += "Not used";
701 l.DrawLatex(0.1,0.6,listOfCuts.Data());
702 listOfCuts = "Maximum cov33: ";
703 if(fMaxCov33Flag) listOfCuts += fMaxCov33;
704 else listOfCuts += "Not used";
705 l.DrawLatex(0.1,0.45,listOfCuts.Data());
706 listOfCuts = "Maximum cov44: ";
707 if(fMaxCov44Flag) listOfCuts += fMaxCov44;
708 else listOfCuts += "Not used";
709 l.DrawLatex(0.1,0.3,listOfCuts.Data());
710 listOfCuts = "Maximum cov55: ";
711 if(fMaxCov55Flag) listOfCuts += fMaxCov55;
712 else listOfCuts += "Not used";
713 l.DrawLatex(0.1,0.15,listOfCuts.Data());
714
715 c->cd(5)->SetFillColor(5);
716 l.DrawLatex(0.3,0.9,"DCA related cuts");
717 listOfCuts = "Maximum sigma to vertex: ";
718 if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
719 else listOfCuts += "Not used";
720 l.DrawLatex(0.1,0.8,listOfCuts.Data());
721 listOfCuts = "Maximum sigma to vertex (TPC): ";
722 if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
723 else listOfCuts += "Not used";
724 l.DrawLatex(0.1,0.72,listOfCuts.Data());
725 listOfCuts = "Maximum DCA in xy: ";
726 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
727 else listOfCuts += "Not used";
728 l.DrawLatex(0.1,0.64,listOfCuts.Data());
729 listOfCuts = "Maximum DCA in z: ";
730 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
731 else listOfCuts += "Not used";
732 l.DrawLatex(0.1,0.56,listOfCuts.Data());
733 listOfCuts = "Maximum DCA in 3D: ";
734 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
735 else listOfCuts += "Not used";
736 l.DrawLatex(0.1,0.48,listOfCuts.Data());
737 listOfCuts = "Maximum DCA in xy (TPC): ";
738 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
739 else listOfCuts += "Not used";
740 l.DrawLatex(0.1,0.4,listOfCuts.Data());
741 listOfCuts = "Maximum DCA in z (TPC): ";
742 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
743 else listOfCuts += "Not used";
744 l.DrawLatex(0.1,0.32,listOfCuts.Data());
745 listOfCuts = "Maximum DCA in 3D (TPC): ";
746 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
747 else listOfCuts += "Not used";
748 l.DrawLatex(0.1,0.24,listOfCuts.Data());
749 listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
750 if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
751 else listOfCuts += "Not used";
752 l.DrawLatex(0.1,0.16,listOfCuts.Data());
753
754 c->cd(6)->SetFillColor(6);
755 l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
756 listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
757 l.DrawLatex(0.1,0.7,listOfCuts.Data());
758 listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
759 l.DrawLatex(0.1,0.5,listOfCuts.Data());
760 listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
761 l.DrawLatex(0.1,0.3,listOfCuts.Data());
762 listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
763 l.DrawLatex(0.1,0.1,listOfCuts.Data());
764
765 return c;
766}
767
768//________________________________________________________________________
769Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
770 //Function that checks if a track is a proton
771 Double_t probability[5];
87a55728 772 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
0ab648ea 773 Long64_t fParticleType = 0;
774
775 //Bayesian approach for the PID
776 if(fProtonPIDMode == kBayesian) {
777 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
778 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
779 if(tpcTrack) {
780 gPt = tpcTrack->Pt();
781 gP = tpcTrack->P();
782 track->GetTPCpid(probability);
783 }
784 }//TPC standalone or Hybrid TPC
785 else if(fProtonAnalysisMode == kGlobal) {
786 gPt = track->Pt();
787 gP = track->P();
788 track->GetESDpid(probability);
789 }//Global tracking
790
791 Double_t rcc = 0.0;
792 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
793 rcc += probability[i]*GetParticleFraction(i,gP);
794 if(rcc != 0.0) {
795 Double_t w[5];
796 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
797 w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
798 fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
799 }
800 if(fParticleType == 4)
801 return kTRUE;
802 }
803 //Ratio of the measured over the theoretical dE/dx a la STAR
804 else if(fProtonPIDMode == kRatio) {
805 Printf("The kRatio mode is not implemented yet!!!");
806 return kFALSE;
807 }
808 //Definition of an N-sigma area around the dE/dx vs P band
87a55728 809 else if(fProtonPIDMode == kSigma1) {
810 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
811 if(tpcTrack) {
812 gPt = tpcTrack->Pt();
813 gP = tpcTrack->P();
814 gEta = tpcTrack->Eta();
815 }
816 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
817 Int_t nbinP = Int_t(gP/0.05) - 6;
818 Double_t tpcSignal = track->GetTPCsignal();
819 Double_t dEdxMean = fdEdxMean[nbinP];
820 Double_t dEdxSigma = fdEdxSigma[nbinP];
821 if((tpcSignal <= dEdxMean + fNSigma*dEdxSigma)&&
822 (tpcSignal <= dEdxMean + fNSigma*dEdxSigma))
823 return kTRUE;
824 }//kSigma1 PID method
825 //Another definition of an N-sigma area around the dE/dx vs P band
826 else if(fProtonPIDMode == kSigma2) {
827 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
828 if(tpcTrack) {
829 gPt = tpcTrack->Pt();
830 gP = tpcTrack->P();
831 gEta = tpcTrack->Eta();
832 }
833 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
834 Int_t nbinP = Int_t(gP/0.05) - 6;
835 Double_t tpcSignal = track->GetTPCsignal();
836 Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
837 Double_t dEdxSigma = fdEdxSigma[nbinP];
838 Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
839 if(nsigma <= fNSigma)
840 return kTRUE;
0ab648ea 841 }
842
843 return kFALSE;
844}
845
87a55728 846//________________________________________________________________________
847void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
848 // This function is used in case the kSigma1 or kSigma2 PID mode is selected
849 // It takes as an argument the name of the ascii file (for the time being)
850 // that is generated as a prior process.
851 // This ascii file has three columns: The min. P value (bins of 50MeV/c)
852 // the mean and the sigma of the dE/dx distributions for protons coming
853 // from a gaussian fit.
854 ifstream in;
855 in.open(filename);
856
857 Double_t gPtMin = 0.0;
858 Int_t iCounter = 0;
859 while(in.good()) {
860 in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
861 if(fDebugMode)
862 Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
863 iCounter += 1;
864 }
865}
866
867//________________________________________________________________________
73aba974 868Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
87a55728 869 // This is the Bethe-Bloch function normalised to 1 at the minimum
870 // We renormalize it based on the MC information
871 // WARNING: To be revised soon!!!
872 // This is just a temporary fix
873 Double_t normalization = 49.2;
874 Double_t bg2=bg*bg;
875 Double_t beta2 = bg2/(1.+ bg2);
876 Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
877 //
878 const Float_t kmeanCorrection =0.1;
879 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
880 bb *= meanCorrection;
881
882 return normalization*bb;
883}
0ab648ea 884
885