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