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