]> 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
libTPCcalib.pkg -
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysisBase.cxx
CommitLineData
0ab648ea 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14/* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
15
16//-----------------------------------------------------------------
17// AliProtonAnalysisBase class
18// This is the class to deal with the proton analysis
19// Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21#include <Riostream.h>
0ab648ea 22#include <TCanvas.h>
23#include <TLatex.h>
73aba974 24#include <TF1.h>
735cc63d 25#include <TList.h>
26#include <TH1F.h>
0ab648ea 27
28#include <AliExternalTrackParam.h>
29#include <AliESDEvent.h>
0ab648ea 30#include <AliPID.h>
0ab648ea 31#include <AliVertexerTracks.h>
10d100d4 32#include <AliESDpid.h>
f81f1b19 33#include <AliTPCPIDResponse.h>
73aba974 34class AliLog;
35class AliESDVertex;
36
37#include "AliProtonAnalysisBase.h"
38
0ab648ea 39ClassImp(AliProtonAnalysisBase)
40
41//____________________________________________________________________//
42AliProtonAnalysisBase::AliProtonAnalysisBase() :
790140ac 43 TObject(), fProtonAnalysisLevel("ESD"), fAnalysisMC(kFALSE),
e56f08ed 44 fTriggerMode(kMB2), kUseOnlineTrigger(kFALSE), kUseOfflineTrigger(kFALSE),
45 fPhysicsSelection(0),
76161fe0 46 fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
0ab648ea 47 fAnalysisEtaMode(kFALSE),
e56f08ed 48 fVxMax(100.), fVyMax(100.), fVzMax(100.), fMinNumOfContributors(0),
0ab648ea 49 fNBinsX(0), fMinX(0), fMaxX(0),
50 fNBinsY(0), fMinY(0), fMaxY(0),
51 fMinTPCClusters(0), fMinITSClusters(0),
52 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
53 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
54 fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
55 fMaxDCAXY(0), fMaxDCAXYTPC(0),
56 fMaxDCAZ(0), fMaxDCAZTPC(0),
57 fMaxDCA3D(0), fMaxDCA3DTPC(0),
87a55728 58 fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
0ab648ea 59 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
60 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
61 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
62 fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
63 fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
64 fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
65 fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
66 fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
67 fMaxConstrainChi2Flag(kFALSE),
68 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
f62e9410 69 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
4787e0ca 70 fPointOnSPDLayersFlag(0), fPointOnSDDLayersFlag(0), fPointOnSSDLayersFlag(0),
0ab648ea 71 fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
72 fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
73 fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
87a55728 74 fMinTPCdEdxPointsFlag(kFALSE),
ff1c9f70 75 fPtDependentDcaXY(0), fPtDependentDcaXYFlag(kFALSE), fNSigmaDCAXY(0.0),
0ab648ea 76 fFunctionProbabilityFlag(kFALSE),
c6909683 77 fNSigma(0), fNRatio(0),
0ab648ea 78 fElectronFunction(0), fMuonFunction(0),
79 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
df201289 80 fDebugMode(kFALSE), fListVertexQA(new TList()) {
0ab648ea 81 //Default constructor
82 for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
c6909683 83 /*for(Int_t i = 0; i < 24; i++) {
87a55728 84 fdEdxMean[i] = 0.0;
85 fdEdxSigma[i] = 0.0;
c6909683 86 }*/
735cc63d 87 fListVertexQA->SetName("fListVertexQA");
88 TH1F *gHistVx = new TH1F("gHistVx",
89 "Vx distribution;V_{x} [cm];Entries",
90 100,-5.,5.);
91 gHistVx->SetFillColor(kRed-2);
92 fListVertexQA->Add(gHistVx);
93 TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
94 "Vx distribution;V_{x} [cm];Entries",
95 100,-5.,5.);
96 fListVertexQA->Add(gHistVxAccepted);
97 TH1F *gHistVy = new TH1F("gHistVy",
98 "Vy distribution;V_{y} [cm];Entries",
99 100,-5.,5.);
100 gHistVy->SetFillColor(kRed-2);
101 fListVertexQA->Add(gHistVy);
102 TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
103 "Vy distribution;V_{y} [cm];Entries",
104 100,-5.,5.);
105 fListVertexQA->Add(gHistVyAccepted);
106 TH1F *gHistVz = new TH1F("gHistVz",
107 "Vz distribution;V_{z} [cm];Entries",
108 100,-25.,25.);
109 gHistVz->SetFillColor(kRed-2);
110 fListVertexQA->Add(gHistVz);
111 TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
112 "Vz distribution;V_{z} [cm];Entries",
113 100,-25.,25.);
114 fListVertexQA->Add(gHistVzAccepted);
115
e56f08ed 116 TH1F *gHistNumberOfContributors = new TH1F("gHistNumberOfContributors",
117 "Number of contributors;N_{contr.};Entries",
118 100,0.,100.);
119 fListVertexQA->Add(gHistNumberOfContributors);
120
735cc63d 121
0ab648ea 122}
123
124//____________________________________________________________________//
125AliProtonAnalysisBase::~AliProtonAnalysisBase() {
126 //Default destructor
127 if(fElectronFunction) delete fElectronFunction;
128 if(fMuonFunction) delete fMuonFunction;
129 if(fPionFunction) delete fPionFunction;
130 if(fKaonFunction) delete fKaonFunction;
131 if(fProtonFunction) delete fProtonFunction;
735cc63d 132 if(fListVertexQA) delete fListVertexQA;
0ab648ea 133}
134
135//____________________________________________________________________//
136Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
137 //Return the a priori probs
138 Double_t partFrac=0;
139 if(fFunctionProbabilityFlag) {
140 if(i == 0) partFrac = fElectronFunction->Eval(p);
141 if(i == 1) partFrac = fMuonFunction->Eval(p);
142 if(i == 2) partFrac = fPionFunction->Eval(p);
143 if(i == 3) partFrac = fKaonFunction->Eval(p);
144 if(i == 4) partFrac = fProtonFunction->Eval(p);
145 }
146 else partFrac = fPartFrac[i];
147
148 return partFrac;
149}
150
151//____________________________________________________________________//
152Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
153 // Checks if the track is outside the analyzed y-Pt phase space
df201289 154 Double_t gP = 0.0, gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
0ab648ea 155 Double_t eta = 0.0;
156
afc2ac17 157 if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid) || (fProtonAnalysisMode == kFullHybrid)) {
0ab648ea 158 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
159 if(!tpcTrack) {
df201289 160 gP = 0.0; gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
0ab648ea 161 }
162 else {
df201289 163 gP = tpcTrack->P();
0ab648ea 164 gPt = tpcTrack->Pt();
165 gPx = tpcTrack->Px();
166 gPy = tpcTrack->Py();
167 gPz = tpcTrack->Pz();
168 eta = tpcTrack->Eta();
169 }
a84984fa 170 }//standalone TPC or Hybrid TPC approaches
0ab648ea 171 else {
df201289 172 gP = track->P();
0ab648ea 173 gPt = track->Pt();
174 gPx = track->Px();
175 gPy = track->Py();
176 gPz = track->Pz();
177 eta = track->Eta();
178 }
179
e56f08ed 180 if((gPt < fMinY) || (gPt > fMaxY)) {
0ab648ea 181 if(fDebugMode)
182 Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
183 return kFALSE;
184 }
e56f08ed 185 if((gP < fMinY) || (gP > fMaxY)) {
186 if(fDebugMode)
187 Printf("IsInPhaseSpace: Track rejected because it has a P value of %lf (accepted interval: %lf - %lf)",gP,fMinY,fMaxY);
188 return kFALSE;
189 }
0ab648ea 190 if(fAnalysisEtaMode) {
191 if((eta < fMinX) || (eta > fMaxX)) {
192 if(fDebugMode)
193 Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
194 return kFALSE;
195 }
196 }
197 else {
198 if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
199 if(fDebugMode)
200 Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
201 return kFALSE;
202 }
203 }
204
205 return kTRUE;
206}
207
208//____________________________________________________________________//
30f87a5b 209Bool_t AliProtonAnalysisBase::IsAccepted(AliESDtrack* track) {
0ab648ea 210 // Checks if the track is excluded from the cuts
0ab648ea 211 Int_t fIdxInt[200];
212 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
213 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
214
215 Float_t chi2PerClusterITS = -1;
216 if (nClustersITS!=0)
217 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
218 Float_t chi2PerClusterTPC = -1;
219 if (nClustersTPC!=0)
220 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
221
222 Double_t extCov[15];
223 track->GetExternalCovariance(extCov);
224
4787e0ca 225 if(fPointOnSPDLayersFlag) {
226 if((!track->HasPointOnITSLayer(0))&&(!track->HasPointOnITSLayer(1))) {
227 if(fDebugMode)
228 Printf("IsAccepted: Track rejected because it doesn't have a point on either SPD layers");
229 return kFALSE;
230 }
231 }
232 if(fPointOnSDDLayersFlag) {
233 if((!track->HasPointOnITSLayer(2))&&(!track->HasPointOnITSLayer(3))) {
234 if(fDebugMode)
235 Printf("IsAccepted: Track rejected because it doesn't have a point on either SDD layers");
236 return kFALSE;
237 }
238 }
239 if(fPointOnSSDLayersFlag) {
240 if((!track->HasPointOnITSLayer(4))&&(!track->HasPointOnITSLayer(5))) {
241 if(fDebugMode)
242 Printf("IsAccepted: Track rejected because it doesn't have a point on either SSD layers");
243 return kFALSE;
244 }
245 }
0ab648ea 246 if(fPointOnITSLayer1Flag) {
247 if(!track->HasPointOnITSLayer(0)) {
248 if(fDebugMode)
249 Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
250 return kFALSE;
251 }
252 }
253 if(fPointOnITSLayer2Flag) {
254 if(!track->HasPointOnITSLayer(1)) {
255 if(fDebugMode)
256 Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
257 return kFALSE;
258 }
259 }
260 if(fPointOnITSLayer3Flag) {
261 if(!track->HasPointOnITSLayer(2)) {
262 if(fDebugMode)
263 Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
264 return kFALSE;
265 }
266 }
267 if(fPointOnITSLayer4Flag) {
268 if(!track->HasPointOnITSLayer(3)) {
269 if(fDebugMode)
270 Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
271 return kFALSE;
272 }
273 }
274 if(fPointOnITSLayer5Flag) {
275 if(!track->HasPointOnITSLayer(4)) {
276 if(fDebugMode)
277 Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
278 return kFALSE;
279 }
280 }
281 if(fPointOnITSLayer6Flag) {
282 if(!track->HasPointOnITSLayer(5)) {
283 if(fDebugMode)
284 Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
285 return kFALSE;
286 }
287 }
288 if(fMinITSClustersFlag) {
289 if(nClustersITS < fMinITSClusters) {
290 if(fDebugMode)
291 Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
292 return kFALSE;
293 }
294 }
295 if(fMaxChi2PerITSClusterFlag) {
296 if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
297 if(fDebugMode)
298 Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
299 return kFALSE;
300 }
301 }
302 if(fMinTPCClustersFlag) {
303 if(nClustersTPC < fMinTPCClusters) {
304 if(fDebugMode)
305 Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
306 return kFALSE;
307 }
308 }
309 if(fMaxChi2PerTPCClusterFlag) {
310 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
311 if(fDebugMode)
312 Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
313 return kFALSE;
314 }
315 }
316 if(fMaxCov11Flag) {
317 if(extCov[0] > fMaxCov11) {
318 if(fDebugMode)
319 Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
320 return kFALSE;
321 }
322 }
323 if(fMaxCov22Flag) {
324 if(extCov[2] > fMaxCov22) {
325 if(fDebugMode)
326 Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
327 return kFALSE;
328 }
329 }
330 if(fMaxCov33Flag) {
331 if(extCov[5] > fMaxCov33) {
332 if(fDebugMode)
333 Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
334 return kFALSE;
335 }
336 }
337 if(fMaxCov44Flag) {
338 if(extCov[9] > fMaxCov44) {
339 if(fDebugMode)
340 Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
341 return kFALSE;
342 }
343 }
344 if(fMaxCov55Flag) {
345 if(extCov[14] > fMaxCov55) {
346 if(fDebugMode)
347 Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
348 return kFALSE;
349 }
350 }
87a55728 351 if(fMinTPCdEdxPointsFlag) {
352 if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
353 if(fDebugMode)
354 Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
355 return kFALSE;
356 }
357 }
0ab648ea 358 if(fITSRefitFlag) {
359 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
360 if(fDebugMode)
361 Printf("IsAccepted: Track rejected because it has no ITS refit flag");
362 return kFALSE;
363 }
364 }
365 if(fTPCRefitFlag) {
366 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
367 if(fDebugMode)
368 Printf("IsAccepted: Track rejected because it has no TPC refit flag");
369 return kFALSE;
370 }
371 }
372 if(fESDpidFlag) {
373 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
374 if(fDebugMode)
375 Printf("IsAccepted: Track rejected because it has no ESD pid flag");
376 return kFALSE;
377 }
378 }
379 if(fTPCpidFlag) {
380 if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
381 if(fDebugMode)
382 Printf("IsAccepted: Track rejected because it has no TPC pid flag");
383 return kFALSE;
384 }
385 }
f62e9410 386 if(fTOFpidFlag) {
387 if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
388 if(fDebugMode)
389 Printf("IsAccepted: Track rejected because it has no TOF pid flag");
390 return kFALSE;
391 }
392 }
0ab648ea 393
394 return kTRUE;
395}
396
ff1c9f70 397//____________________________________________________________________//
398void AliProtonAnalysisBase::SetPtDependentDCAxy(Int_t nSigma,
399 Double_t p0,
400 Double_t p1,
401 Double_t p2) {
402 //Pt dependent dca cut in xy
403 fNSigmaDCAXY = nSigma;
404 fPtDependentDcaXY = new TF1("fPtDependentDcaXY","[0]+[1]/x^[2]",0.1,10.1);
405 fPtDependentDcaXY->SetParameter(0,p0);
406 fPtDependentDcaXY->SetParameter(1,p1);
407 fPtDependentDcaXY->SetParameter(2,p2);
408 fPtDependentDcaXYFlag = kTRUE;
409}
410
e56f08ed 411//____________________________________________________________________//
412Bool_t AliProtonAnalysisBase::IsPrimary(AliESDEvent *esd,
413 const AliESDVertex *vertex,
414 AliESDtrack* track) {
415 // Checks if the track is a primary-like candidate
842c574e 416 const Double_t kMicrometer2Centimeter = 0.0001;
e56f08ed 417 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
418 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
419 Double_t dca3D = 0.0;
afc2ac17 420 Float_t dcaXY = 0.0, dcaZ = 0.0;
421
e56f08ed 422 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
423 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
424 if(!tpcTrack) {
425 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
426 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
427 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
428 }
429 else {
430 gPt = tpcTrack->Pt();
431 gPx = tpcTrack->Px();
432 gPy = tpcTrack->Py();
433 gPz = tpcTrack->Pz();
434 tpcTrack->PropagateToDCA(vertex,
435 esd->GetMagneticField(),
436 100.,dca,cov);
437 }
438 }//standalone TPC or hybrid TPC approaches
afc2ac17 439 if(fProtonAnalysisMode == kFullHybrid) {
440 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
441 AliExternalTrackParam cParam;
442 if(!tpcTrack) {
443 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
444 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
445 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
446 }
447 else {
448 gPt = tpcTrack->Pt();
449 gPx = tpcTrack->Px();
450 gPy = tpcTrack->Py();
451 gPz = tpcTrack->Pz();
452 track->RelateToVertex(vertex,
453 esd->GetMagneticField(),
454 100.,&cParam);
455 track->GetImpactParameters(dcaXY,dcaZ);
456 dca[0] = dcaXY; dca[1] = dcaZ;
457 }
458 }//standalone TPC or hybrid TPC approaches
e56f08ed 459 else {
460 gPt = track->Pt();
461 gPx = track->Px();
462 gPy = track->Py();
463 gPz = track->Pz();
464 track->PropagateToDCA(vertex,
465 esd->GetMagneticField(),
466 100.,dca,cov);
467 }
468 dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
469 TMath::Power(dca[1],2));
470
e56f08ed 471 if(fMaxSigmaToVertexFlag) {
472 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
473 if(fDebugMode)
474 Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
475 return kFALSE;
476 }
477 }
478 if(fMaxSigmaToVertexTPCFlag) {
479 if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
480 if(fDebugMode)
481 Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
482 return kFALSE;
483 }
484 }
485 if(fMaxDCAXYFlag) {
486 if(TMath::Abs(dca[0]) > fMaxDCAXY) {
487 if(fDebugMode)
488 Printf("IsPrimary: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
489 return kFALSE;
490 }
491 }
492 if(fMaxDCAXYTPCFlag) {
493 if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
494 if(fDebugMode)
495 Printf("IsPrimary: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
496 return kFALSE;
497 }
498 }
499 if(fMaxDCAZFlag) {
500 if(TMath::Abs(dca[1]) > fMaxDCAZ) {
501 if(fDebugMode)
502 Printf("IsPrimary: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
503 return kFALSE;
504 }
505 }
506 if(fMaxDCAZTPCFlag) {
507 if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
508 if(fDebugMode)
509 Printf("IsPrimary: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
510 return kFALSE;
511 }
512 }
513 if(fMaxDCA3DFlag) {
514 if(TMath::Abs(dca3D) > fMaxDCA3D) {
515 if(fDebugMode)
516 Printf("IsPrimary: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
517 return kFALSE;
518 }
519 }
520 if(fMaxDCA3DTPCFlag) {
521 if(TMath::Abs(dca3D) > fMaxDCA3DTPC) {
522 if(fDebugMode)
523 Printf("IsPrimary: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
524 return kFALSE;
525 }
526 }
527 if(fMaxConstrainChi2Flag) {
528 if(track->GetConstrainedChi2() > 0)
529 if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) {
530 if(fDebugMode)
531 Printf("IsPrimary: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
532 return kFALSE;
533 }
534 }
ff1c9f70 535 if(fPtDependentDcaXYFlag) {
842c574e 536 if(TMath::Abs(dca[0]) > kMicrometer2Centimeter*fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt)) {
ff1c9f70 537 if(fDebugMode)
538 Printf("IsPrimary: Track rejected because it has a value of the dca(xy) higher than the %d sigma pt dependent cut: %lf (max. requested: %lf)",TMath::Abs(dca[0]),fNSigmaDCAXY,fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt));
539 return kFALSE;
540 }
541 }
e56f08ed 542
543 return kTRUE;
544}
545
0ab648ea 546//____________________________________________________________________//
547Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
548 // Calculates the number of sigma to the vertex.
0ab648ea 549 Float_t b[2];
550 Float_t bRes[2];
551 Float_t bCov[3];
afc2ac17 552 if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid)&&(fProtonAnalysisMode != kFullHybrid))
0ab648ea 553 esdTrack->GetImpactParametersTPC(b,bCov);
554 else
555 esdTrack->GetImpactParameters(b,bCov);
556
557 if (bCov[0]<=0 || bCov[2]<=0) {
558 //AliDebug(1, "Estimated b resolution lower or equal zero!");
559 bCov[0]=0; bCov[2]=0;
560 }
561 bRes[0] = TMath::Sqrt(bCov[0]);
562 bRes[1] = TMath::Sqrt(bCov[2]);
563
564 if (bRes[0] == 0 || bRes[1] ==0) return -1;
565
566 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
567
568 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
569
570 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
571
572 return d;
573}
574
575//____________________________________________________________________//
576Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
577 Double_t gPy,
578 Double_t gPz) const {
579 //returns the rapidity of the proton - to be removed
580 Double_t fMass = 9.38270000000000048e-01;
581
582 Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
583 TMath::Power(gPy,2) +
584 TMath::Power(gPz,2));
585 Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
586 Double_t y = -999;
587 if(energy != gPz)
588 y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
589
590 return y;
591}
592
593//________________________________________________________________________
594const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
595 AnalysisMode mode,
596 Double_t gVxMax,
597 Double_t gVyMax,
598 Double_t gVzMax) {
599 // Get the vertex from the ESD and returns it if the vertex is valid
600 // Second argument decides which vertex is used (this selects
601 // also the quality criteria that are applied)
602 const AliESDVertex* vertex = 0;
afc2ac17 603 if((mode == kHybrid)||(mode == kFullHybrid))
0ab648ea 604 vertex = esd->GetPrimaryVertexSPD();
605 else if(mode == kTPC){
606 Double_t kBz = esd->GetMagneticField();
607 AliVertexerTracks vertexer(kBz);
608 vertexer.SetTPCMode();
609 AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
610 esd->SetPrimaryVertexTPC(vTPC);
611 for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
612 AliESDtrack *t = esd->GetTrack(i);
613 t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
614 }
615 delete vTPC;
616 vertex = esd->GetPrimaryVertexTPC();
617 }
618 else if(mode == kGlobal)
619 vertex = esd->GetPrimaryVertex();
620 else
621 Printf("GetVertex: ERROR: Invalid second argument %d", mode);
622
623 if(!vertex) {
624 if(fDebugMode)
625 Printf("GetVertex: Event rejected because there is no valid vertex object");
626 return 0;
627 }
628
629 // check Ncontributors
630 if(vertex->GetNContributors() <= 0) {
631 if(fDebugMode)
632 Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
633 return 0;
634 }
635
636 // check resolution
637 Double_t zRes = vertex->GetZRes();
638 if(zRes == 0) {
639 if(fDebugMode)
640 Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
641 return 0;
642 }
735cc63d 643 ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
644 ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
645 ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
e56f08ed 646
0ab648ea 647 //check position
648 if(TMath::Abs(vertex->GetXv()) > gVxMax) {
649 if(fDebugMode)
650 Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
651 return 0;
652 }
653 if(TMath::Abs(vertex->GetYv()) > gVyMax) {
654 if(fDebugMode)
655 Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
656 return 0;
657 }
658 if(TMath::Abs(vertex->GetZv()) > gVzMax) {
659 if(fDebugMode)
660 Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
661 return 0;
662 }
735cc63d 663 ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
664 ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
665 ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
e56f08ed 666 ((TH1F *)(fListVertexQA->At(6)))->Fill(vertex->GetNContributors());
667
668 //check number of contributors
669 if(fMinNumOfContributors > 0) {
670 if(fMinNumOfContributors > vertex->GetNContributors()) {
671 if(fDebugMode)
672 Printf("GetVertex: Event rejected because it has %d number of contributors (requested minimum: %d)",vertex->GetNContributors(),fMinNumOfContributors);
673
674 return 0;
675 }
676 }
0ab648ea 677
678 return vertex;
679}
680
681//________________________________________________________________________
682Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
683 TriggerMode trigger) {
684 // check if the event was triggered
685 ULong64_t triggerMask = esd->GetTriggerMask();
42270c4c 686 TString firedTriggerClass = esd->GetFiredTriggerClasses();
0ab648ea 687
42270c4c 688 if(fAnalysisMC) {
689 // definitions from p-p.cfg
690 ULong64_t spdFO = (1 << 14);
691 ULong64_t v0left = (1 << 11);
692 ULong64_t v0right = (1 << 12);
693
694 switch (trigger) {
695 case kMB1: {
696 if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
697 return kTRUE;
698 break;
699 }
700 case kMB2: {
701 if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
702 return kTRUE;
703 break;
704 }
705 case kSPDFASTOR: {
706 if (triggerMask & spdFO)
707 return kTRUE;
708 break;
709 }
710 }//switch
0ab648ea 711 }
42270c4c 712 else {
4de4661f 713 if(kUseOnlineTrigger) {
714 if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
715 return kTRUE;
716 }
717 else if(!kUseOnlineTrigger)
0ab648ea 718 return kTRUE;
0ab648ea 719 }
0ab648ea 720
721 return kFALSE;
722}
723
724//________________________________________________________________________
725TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
726 // return the list of cuts and their cut values for reference
727 TLatex l;
728 l.SetTextAlign(12);
729 l.SetTextSize(0.04);
730
731 TString listOfCuts;
732
733 TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
734 c->SetFillColor(10); c->SetHighLightColor(41);
735 c->Divide(3,2);
736
737 c->cd(1)->SetFillColor(10);
738 l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
739
740 listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
741 l.DrawLatex(0.1,0.82,listOfCuts.Data());
742 listOfCuts = "Analysis mode: ";
743 if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
744 if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
afc2ac17 745 if(fProtonAnalysisMode == kFullHybrid) listOfCuts += "Full Hybrid TPC";
0ab648ea 746 if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
747 l.DrawLatex(0.1,0.74,listOfCuts.Data());
748 listOfCuts = "Trigger mode: ";
749 if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
750 if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
751 if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
752 l.DrawLatex(0.1,0.66,listOfCuts.Data());
753 listOfCuts = "PID mode: ";
754 if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
c6909683 755 if(fProtonPIDMode == kRatio) {
756 listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.}) > ";
757 listOfCuts += fNRatio;
87a55728 758 }
c6909683 759 if(fProtonPIDMode == kSigma) {
760 listOfCuts += "N_{#sigma} area: "; listOfCuts += fNSigma;
87a55728 761 listOfCuts += " #sigma";
762 }
c6909683 763 //if(fProtonPIDMode == kSigma2) {
764 //listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
765 //listOfCuts += " #sigma";
766 //}
0ab648ea 767 l.DrawLatex(0.1,0.58,listOfCuts.Data());
768 listOfCuts = "Accepted vertex diamond: ";
c6909683 769 l.DrawLatex(0.1,0.52,listOfCuts.Data());
0ab648ea 770 listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
c6909683 771 l.DrawLatex(0.6,0.52,listOfCuts.Data());
0ab648ea 772 listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
c6909683 773 l.DrawLatex(0.6,0.45,listOfCuts.Data());
0ab648ea 774 listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
c6909683 775 l.DrawLatex(0.6,0.38,listOfCuts.Data());
776 listOfCuts = "N_{contributors} > "; listOfCuts += fMinNumOfContributors;
777 l.DrawLatex(0.6,0.31,listOfCuts.Data());
0ab648ea 778 listOfCuts = "Phase space: ";
779 l.DrawLatex(0.1,0.2,listOfCuts.Data());
780 if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
781 else listOfCuts = "|y| < ";
782 listOfCuts += TMath::Abs(fMinX);
783 l.DrawLatex(0.35,0.2,listOfCuts.Data());
784 listOfCuts = "N_{bins} = ";
785 listOfCuts += fNBinsX; listOfCuts += " (binning: ";
786 listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
787 l.DrawLatex(0.35,0.15,listOfCuts.Data());
788 listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
789 listOfCuts += fMaxY; listOfCuts += "GeV/c";
790 l.DrawLatex(0.35,0.1,listOfCuts.Data());
791 listOfCuts = "N_{bins} = ";
792 listOfCuts += fNBinsY; listOfCuts += " (binning: ";
793 listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
794 l.DrawLatex(0.35,0.05,listOfCuts.Data());
795
796 c->cd(2)->SetFillColor(2);
4787e0ca 797 l.DrawLatex(0.3,0.95,"ITS related cuts");
798 listOfCuts = "Request a cluster on either of the SPD layers: ";
799 listOfCuts += fPointOnSPDLayersFlag;
800 l.DrawLatex(0.1,0.9,listOfCuts.Data());
801 listOfCuts = "Request a cluster on either of the SDD layers: ";
802 listOfCuts += fPointOnSDDLayersFlag;
803 l.DrawLatex(0.1,0.83,listOfCuts.Data());
804 listOfCuts = "Request a cluster on either of the SSD layers: ";
805 listOfCuts += fPointOnSSDLayersFlag;
806 l.DrawLatex(0.1,0.76,listOfCuts.Data());
807
0ab648ea 808 listOfCuts = "Request a cluster on SPD1: ";
809 listOfCuts += fPointOnITSLayer1Flag;
4787e0ca 810 l.DrawLatex(0.1,0.69,listOfCuts.Data());
0ab648ea 811 listOfCuts = "Request a cluster on SPD2: ";
812 listOfCuts += fPointOnITSLayer2Flag;
4787e0ca 813 l.DrawLatex(0.1,0.62,listOfCuts.Data());
0ab648ea 814 listOfCuts = "Request a cluster on SDD1: ";
815 listOfCuts += fPointOnITSLayer3Flag;
4787e0ca 816 l.DrawLatex(0.1,0.55,listOfCuts.Data());
0ab648ea 817 listOfCuts = "Request a cluster on SDD2: ";
818 listOfCuts += fPointOnITSLayer4Flag;
4787e0ca 819 l.DrawLatex(0.1,0.48,listOfCuts.Data());
0ab648ea 820 listOfCuts = "Request a cluster on SSD1: ";
821 listOfCuts += fPointOnITSLayer5Flag;
4787e0ca 822 l.DrawLatex(0.1,0.41,listOfCuts.Data());
0ab648ea 823 listOfCuts = "Request a cluster on SSD2: ";
824 listOfCuts += fPointOnITSLayer6Flag;
4787e0ca 825 l.DrawLatex(0.1,0.34,listOfCuts.Data());
0ab648ea 826 listOfCuts = "Minimum number of ITS clusters: ";
827 if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
828 else listOfCuts += "Not used";
4787e0ca 829 l.DrawLatex(0.1,0.27,listOfCuts.Data());
0ab648ea 830 listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
831 if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
832 else listOfCuts += "Not used";
4787e0ca 833 l.DrawLatex(0.1,0.2,listOfCuts.Data());
0ab648ea 834
835 c->cd(3)->SetFillColor(3);
836 l.DrawLatex(0.3,0.9,"TPC related cuts");
837 listOfCuts = "Minimum number of TPC clusters: ";
838 if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
839 else listOfCuts += "Not used";
840 l.DrawLatex(0.1,0.7,listOfCuts.Data());
841 listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
842 if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
843 else listOfCuts += "Not used";
844 l.DrawLatex(0.1,0.5,listOfCuts.Data());
57e749bb 845 listOfCuts = "Minimum number of TPC points for the dE/dx: ";
846 if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
847 else listOfCuts += "Not used";
848 l.DrawLatex(0.1,0.3,listOfCuts.Data());
0ab648ea 849
850 c->cd(4)->SetFillColor(4);
851 l.DrawLatex(0.3,0.9,"Tracking related cuts");
852 listOfCuts = "Maximum cov11: ";
853 if(fMaxCov11Flag) listOfCuts += fMaxCov11;
854 else listOfCuts += "Not used";
855 l.DrawLatex(0.1,0.75,listOfCuts.Data());
856 listOfCuts = "Maximum cov22: ";
857 if(fMaxCov22Flag) listOfCuts += fMaxCov22;
858 else listOfCuts += "Not used";
859 l.DrawLatex(0.1,0.6,listOfCuts.Data());
860 listOfCuts = "Maximum cov33: ";
861 if(fMaxCov33Flag) listOfCuts += fMaxCov33;
862 else listOfCuts += "Not used";
863 l.DrawLatex(0.1,0.45,listOfCuts.Data());
864 listOfCuts = "Maximum cov44: ";
865 if(fMaxCov44Flag) listOfCuts += fMaxCov44;
866 else listOfCuts += "Not used";
867 l.DrawLatex(0.1,0.3,listOfCuts.Data());
868 listOfCuts = "Maximum cov55: ";
869 if(fMaxCov55Flag) listOfCuts += fMaxCov55;
870 else listOfCuts += "Not used";
871 l.DrawLatex(0.1,0.15,listOfCuts.Data());
872
873 c->cd(5)->SetFillColor(5);
874 l.DrawLatex(0.3,0.9,"DCA related cuts");
875 listOfCuts = "Maximum sigma to vertex: ";
876 if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
877 else listOfCuts += "Not used";
878 l.DrawLatex(0.1,0.8,listOfCuts.Data());
879 listOfCuts = "Maximum sigma to vertex (TPC): ";
880 if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
881 else listOfCuts += "Not used";
882 l.DrawLatex(0.1,0.72,listOfCuts.Data());
883 listOfCuts = "Maximum DCA in xy: ";
884 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
885 else listOfCuts += "Not used";
886 l.DrawLatex(0.1,0.64,listOfCuts.Data());
887 listOfCuts = "Maximum DCA in z: ";
888 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
889 else listOfCuts += "Not used";
890 l.DrawLatex(0.1,0.56,listOfCuts.Data());
891 listOfCuts = "Maximum DCA in 3D: ";
892 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
893 else listOfCuts += "Not used";
894 l.DrawLatex(0.1,0.48,listOfCuts.Data());
895 listOfCuts = "Maximum DCA in xy (TPC): ";
896 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
897 else listOfCuts += "Not used";
898 l.DrawLatex(0.1,0.4,listOfCuts.Data());
899 listOfCuts = "Maximum DCA in z (TPC): ";
900 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
901 else listOfCuts += "Not used";
902 l.DrawLatex(0.1,0.32,listOfCuts.Data());
903 listOfCuts = "Maximum DCA in 3D (TPC): ";
904 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
905 else listOfCuts += "Not used";
906 l.DrawLatex(0.1,0.24,listOfCuts.Data());
907 listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
908 if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
909 else listOfCuts += "Not used";
910 l.DrawLatex(0.1,0.16,listOfCuts.Data());
911
912 c->cd(6)->SetFillColor(6);
913 l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
914 listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
915 l.DrawLatex(0.1,0.7,listOfCuts.Data());
916 listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
917 l.DrawLatex(0.1,0.5,listOfCuts.Data());
918 listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
919 l.DrawLatex(0.1,0.3,listOfCuts.Data());
920 listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
921 l.DrawLatex(0.1,0.1,listOfCuts.Data());
922
923 return c;
924}
925
926//________________________________________________________________________
927Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
928 //Function that checks if a track is a proton
929 Double_t probability[5];
87a55728 930 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
0ab648ea 931 Long64_t fParticleType = 0;
932
933 //Bayesian approach for the PID
934 if(fProtonPIDMode == kBayesian) {
afc2ac17 935 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)||(fProtonAnalysisMode == kFullHybrid)) {
0ab648ea 936 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
937 if(tpcTrack) {
938 gPt = tpcTrack->Pt();
939 gP = tpcTrack->P();
940 track->GetTPCpid(probability);
941 }
942 }//TPC standalone or Hybrid TPC
943 else if(fProtonAnalysisMode == kGlobal) {
944 gPt = track->Pt();
945 gP = track->P();
946 track->GetESDpid(probability);
947 }//Global tracking
948
949 Double_t rcc = 0.0;
950 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
951 rcc += probability[i]*GetParticleFraction(i,gP);
952 if(rcc != 0.0) {
953 Double_t w[5];
954 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
955 w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
956 fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
957 }
958 if(fParticleType == 4)
959 return kTRUE;
960 }
961 //Ratio of the measured over the theoretical dE/dx a la STAR
962 else if(fProtonPIDMode == kRatio) {
c6909683 963 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
964 if(tpcTrack) {
965 gPt = tpcTrack->Pt();
b8cca027 966 gP = track->GetInnerParam()->P();
c6909683 967 gEta = tpcTrack->Eta();
968 }
969 Double_t fAlephParameters[5];
970 if(fAnalysisMC) {
971 fAlephParameters[0] = 2.15898e+00/50.;
972 fAlephParameters[1] = 1.75295e+01;
973 fAlephParameters[2] = 3.40030e-09;
974 fAlephParameters[3] = 1.96178e+00;
975 fAlephParameters[4] = 3.91720e+00;
976 }
977 else {
978 fAlephParameters[0] = 0.0283086;
979 fAlephParameters[1] = 2.63394e+01;
980 fAlephParameters[2] = 5.04114e-11;
981 fAlephParameters[3] = 2.12543e+00;
982 fAlephParameters[4] = 4.88663e+00;
983 }
984
985 AliESDpid *fESDpid = new AliESDpid();
986 AliTPCPIDResponse tpcResponse = fESDpid->GetTPCResponse();
987 tpcResponse.SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
9a3204e2 988
989 Double_t normalizeddEdx = -10.;
990 if((track->GetTPCsignal() > 0.0) && (tpcResponse.GetExpectedSignal(gP,AliPID::kProton) > 0.0))
991 TMath::Log(track->GetTPCsignal()/tpcResponse.GetExpectedSignal(gP,AliPID::kProton));
c6909683 992
993 if(normalizeddEdx >= fNRatio)
994 return kTRUE;
995 }//kRatio PID mode
0ab648ea 996 //Definition of an N-sigma area around the dE/dx vs P band
c6909683 997 else if(fProtonPIDMode == kSigma) {
2c080079 998 Double_t fAlephParameters[5];
0bb8f6b4 999 if(fAnalysisMC) {
e56f08ed 1000 fAlephParameters[0] = 2.15898e+00/50.;
1001 fAlephParameters[1] = 1.75295e+01;
1002 fAlephParameters[2] = 3.40030e-09;
1003 fAlephParameters[3] = 1.96178e+00;
1004 fAlephParameters[4] = 3.91720e+00;
0bb8f6b4 1005 }
1006 else {
e56f08ed 1007 fAlephParameters[0] = 0.0283086;
1008 fAlephParameters[1] = 2.63394e+01;
1009 fAlephParameters[2] = 5.04114e-11;
1010 fAlephParameters[3] = 2.12543e+00;
1011 fAlephParameters[4] = 4.88663e+00;
0bb8f6b4 1012 }
2c080079 1013
30f87a5b 1014 Double_t nsigma = 100.0;
10d100d4 1015 AliESDpid *fESDpid = new AliESDpid();
e56f08ed 1016 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
2c080079 1017
30f87a5b 1018 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
1019 if(tpcTrack)
1020 nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton));
1021
2c080079 1022 if(nsigma <= fNSigma)
87a55728 1023 return kTRUE;
c6909683 1024 }//kSigma PID method
87a55728 1025 //Another definition of an N-sigma area around the dE/dx vs P band
c6909683 1026 /*else if(fProtonPIDMode == kSigma2) {
87a55728 1027 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
1028 if(tpcTrack) {
1029 gPt = tpcTrack->Pt();
1030 gP = tpcTrack->P();
1031 gEta = tpcTrack->Eta();
1032 }
f81f1b19 1033 Double_t fAlephParameters[5];
1034 if(fAnalysisMC) {
1035 fAlephParameters[0] = 2.15898e+00/50.;
1036 fAlephParameters[1] = 1.75295e+01;
1037 fAlephParameters[2] = 3.40030e-09;
1038 fAlephParameters[3] = 1.96178e+00;
1039 fAlephParameters[4] = 3.91720e+00;
1040 }
1041 else {
1042 fAlephParameters[0] = 0.0283086;
1043 fAlephParameters[1] = 2.63394e+01;
1044 fAlephParameters[2] = 5.04114e-11;
1045 fAlephParameters[3] = 2.12543e+00;
1046 fAlephParameters[4] = 4.88663e+00;
1047 }
1048
1049 AliESDpid *fESDpid = new AliESDpid();
1050 AliTPCPIDResponse tpcResponse = fESDpid->GetTPCResponse();
1051 tpcResponse.SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
1052 Double_t normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse.GetExpectedSignal(gP,AliPID::kProton));
1053
1054 if(normalizeddEdx >= -0.15)
87a55728 1055 return kTRUE;
c6909683 1056 }*/
0ab648ea 1057
1058 return kFALSE;
1059}
1060
1061
1062