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