]> 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
Addition of misaligner classes and small changes
[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());
57e749bb 691 listOfCuts = "Minimum number of TPC points for the dE/dx: ";
692 if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
693 else listOfCuts += "Not used";
694 l.DrawLatex(0.1,0.3,listOfCuts.Data());
0ab648ea 695
696 c->cd(4)->SetFillColor(4);
697 l.DrawLatex(0.3,0.9,"Tracking related cuts");
698 listOfCuts = "Maximum cov11: ";
699 if(fMaxCov11Flag) listOfCuts += fMaxCov11;
700 else listOfCuts += "Not used";
701 l.DrawLatex(0.1,0.75,listOfCuts.Data());
702 listOfCuts = "Maximum cov22: ";
703 if(fMaxCov22Flag) listOfCuts += fMaxCov22;
704 else listOfCuts += "Not used";
705 l.DrawLatex(0.1,0.6,listOfCuts.Data());
706 listOfCuts = "Maximum cov33: ";
707 if(fMaxCov33Flag) listOfCuts += fMaxCov33;
708 else listOfCuts += "Not used";
709 l.DrawLatex(0.1,0.45,listOfCuts.Data());
710 listOfCuts = "Maximum cov44: ";
711 if(fMaxCov44Flag) listOfCuts += fMaxCov44;
712 else listOfCuts += "Not used";
713 l.DrawLatex(0.1,0.3,listOfCuts.Data());
714 listOfCuts = "Maximum cov55: ";
715 if(fMaxCov55Flag) listOfCuts += fMaxCov55;
716 else listOfCuts += "Not used";
717 l.DrawLatex(0.1,0.15,listOfCuts.Data());
718
719 c->cd(5)->SetFillColor(5);
720 l.DrawLatex(0.3,0.9,"DCA related cuts");
721 listOfCuts = "Maximum sigma to vertex: ";
722 if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
723 else listOfCuts += "Not used";
724 l.DrawLatex(0.1,0.8,listOfCuts.Data());
725 listOfCuts = "Maximum sigma to vertex (TPC): ";
726 if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
727 else listOfCuts += "Not used";
728 l.DrawLatex(0.1,0.72,listOfCuts.Data());
729 listOfCuts = "Maximum DCA in xy: ";
730 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
731 else listOfCuts += "Not used";
732 l.DrawLatex(0.1,0.64,listOfCuts.Data());
733 listOfCuts = "Maximum DCA in z: ";
734 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
735 else listOfCuts += "Not used";
736 l.DrawLatex(0.1,0.56,listOfCuts.Data());
737 listOfCuts = "Maximum DCA in 3D: ";
738 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
739 else listOfCuts += "Not used";
740 l.DrawLatex(0.1,0.48,listOfCuts.Data());
741 listOfCuts = "Maximum DCA in xy (TPC): ";
742 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
743 else listOfCuts += "Not used";
744 l.DrawLatex(0.1,0.4,listOfCuts.Data());
745 listOfCuts = "Maximum DCA in z (TPC): ";
746 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
747 else listOfCuts += "Not used";
748 l.DrawLatex(0.1,0.32,listOfCuts.Data());
749 listOfCuts = "Maximum DCA in 3D (TPC): ";
750 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
751 else listOfCuts += "Not used";
752 l.DrawLatex(0.1,0.24,listOfCuts.Data());
753 listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
754 if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
755 else listOfCuts += "Not used";
756 l.DrawLatex(0.1,0.16,listOfCuts.Data());
757
758 c->cd(6)->SetFillColor(6);
759 l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
760 listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
761 l.DrawLatex(0.1,0.7,listOfCuts.Data());
762 listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
763 l.DrawLatex(0.1,0.5,listOfCuts.Data());
764 listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
765 l.DrawLatex(0.1,0.3,listOfCuts.Data());
766 listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
767 l.DrawLatex(0.1,0.1,listOfCuts.Data());
768
769 return c;
770}
771
772//________________________________________________________________________
773Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
774 //Function that checks if a track is a proton
775 Double_t probability[5];
87a55728 776 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
0ab648ea 777 Long64_t fParticleType = 0;
778
779 //Bayesian approach for the PID
780 if(fProtonPIDMode == kBayesian) {
781 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
782 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
783 if(tpcTrack) {
784 gPt = tpcTrack->Pt();
785 gP = tpcTrack->P();
786 track->GetTPCpid(probability);
787 }
788 }//TPC standalone or Hybrid TPC
789 else if(fProtonAnalysisMode == kGlobal) {
790 gPt = track->Pt();
791 gP = track->P();
792 track->GetESDpid(probability);
793 }//Global tracking
794
795 Double_t rcc = 0.0;
796 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
797 rcc += probability[i]*GetParticleFraction(i,gP);
798 if(rcc != 0.0) {
799 Double_t w[5];
800 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
801 w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
802 fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
803 }
804 if(fParticleType == 4)
805 return kTRUE;
806 }
807 //Ratio of the measured over the theoretical dE/dx a la STAR
808 else if(fProtonPIDMode == kRatio) {
809 Printf("The kRatio mode is not implemented yet!!!");
810 return kFALSE;
811 }
812 //Definition of an N-sigma area around the dE/dx vs P band
87a55728 813 else if(fProtonPIDMode == kSigma1) {
814 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
815 if(tpcTrack) {
816 gPt = tpcTrack->Pt();
817 gP = tpcTrack->P();
818 gEta = tpcTrack->Eta();
819 }
820 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
821 Int_t nbinP = Int_t(gP/0.05) - 6;
822 Double_t tpcSignal = track->GetTPCsignal();
823 Double_t dEdxMean = fdEdxMean[nbinP];
824 Double_t dEdxSigma = fdEdxSigma[nbinP];
825 if((tpcSignal <= dEdxMean + fNSigma*dEdxSigma)&&
826 (tpcSignal <= dEdxMean + fNSigma*dEdxSigma))
827 return kTRUE;
828 }//kSigma1 PID method
829 //Another definition of an N-sigma area around the dE/dx vs P band
830 else if(fProtonPIDMode == kSigma2) {
831 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
832 if(tpcTrack) {
833 gPt = tpcTrack->Pt();
834 gP = tpcTrack->P();
835 gEta = tpcTrack->Eta();
836 }
837 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
838 Int_t nbinP = Int_t(gP/0.05) - 6;
839 Double_t tpcSignal = track->GetTPCsignal();
840 Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
841 Double_t dEdxSigma = fdEdxSigma[nbinP];
842 Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
843 if(nsigma <= fNSigma)
844 return kTRUE;
0ab648ea 845 }
846
847 return kFALSE;
848}
849
87a55728 850//________________________________________________________________________
851void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
852 // This function is used in case the kSigma1 or kSigma2 PID mode is selected
853 // It takes as an argument the name of the ascii file (for the time being)
854 // that is generated as a prior process.
855 // This ascii file has three columns: The min. P value (bins of 50MeV/c)
856 // the mean and the sigma of the dE/dx distributions for protons coming
857 // from a gaussian fit.
858 ifstream in;
859 in.open(filename);
860
861 Double_t gPtMin = 0.0;
862 Int_t iCounter = 0;
863 while(in.good()) {
864 in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
865 if(fDebugMode)
866 Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
867 iCounter += 1;
868 }
869}
870
871//________________________________________________________________________
73aba974 872Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
87a55728 873 // This is the Bethe-Bloch function normalised to 1 at the minimum
874 // We renormalize it based on the MC information
875 // WARNING: To be revised soon!!!
876 // This is just a temporary fix
877 Double_t normalization = 49.2;
878 Double_t bg2=bg*bg;
879 Double_t beta2 = bg2/(1.+ bg2);
880 Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
881 //
882 const Float_t kmeanCorrection =0.1;
883 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
884 bb *= meanCorrection;
885
886 return normalization*bb;
887}
0ab648ea 888
889