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