1 /*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
19 // Class for impact parameter (DCA) of charged particles
20 // Study resolution and pull: prepare for beauty study
23 // Hongyan Yang <hongyan@physi.uni-heidelberg.de>
24 // Carlo Bombonati <carlo.bombonati@cern.ch>
30 #include <TParticle.h>
31 #include "AliMCParticle.h"
32 #include "AliESDtrack.h"
33 #include "AliESDEvent.h"
34 #include "AliMCEvent.h"
35 #include "AliMCVertex.h"
37 #include "AliKFParticle.h"
38 #include "AliKFVertex.h"
40 #include "AliESDVertex.h"
44 #include "AliHFEdca.h"
49 //________________________________________________________________________
50 const Char_t* AliHFEdca::fgkParticles[12] = {
52 "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
53 "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
54 "allNegative", "allPositive"
57 const Int_t AliHFEdca::fgkPdgParticle[10] = {
58 // 11, 13, -211, -233, -2122,
59 // -11, -13, 211, 233, 2122};
60 kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
61 -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
63 //________________________________________________________________________
64 const Int_t AliHFEdca::fgkColorPart[12] = {
65 // colors assigned to particles
66 kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
67 kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
71 //________________________________________________________________________
72 const Float_t AliHFEdca::fgkPtIntv[51] = {
74 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55,
75 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50,
76 1.70, 1.90, 2.10, 2.30, 2.50, 2.70, 2.90, 3.10, 3.30, 3.50,
77 3.80, 4.10, 4.40, 4.70, 5.00, 5.30, 5.60, 5.90, 6.20, 6.50,
78 7.00, 7.50, 8.00, 9.00, 10.0, 11.0, 12.0, 13.0, 15.0, 18.0,
81 //________________________________________________________________________
82 const Char_t *AliHFEdca::fgkDcaVar[2] = {
85 //________________________________________________________________________
86 const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
87 ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
89 //________________________________________________________________________
90 const Char_t *AliHFEdca::fgkVertexVar[3] = {
91 "VertexX", "VertexY", "VertexZ"};
93 //________________________________________________________________________
94 const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
95 ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
97 //________________________________________________________________________
98 const Char_t *AliHFEdca::fgkResDcaVar[2] = {
99 "deltaDcaXY", "deltaDcaZ"};
101 //________________________________________________________________________
102 const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
103 ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
105 //________________________________________________________________________
106 const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
107 "pullDcaXY", "pullDcaZ"
110 //________________________________________________________________________
111 const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
112 ";residual dca_{xy}/(error dca_{xy});counts",
113 ";residual dca_{z}/(error dca_{z});counts"
116 //________________________________________________________________________
117 const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
118 ";dca_{xy}^{data}/error dca_{xy};counts",
119 ";dca_{z}^{data}/error dca_{z};counts"
122 //________________________________________________________________________
123 AliHFEdca::AliHFEdca():
127 // default constructor
129 for(Int_t j=0; j<kNParticles; j++){
130 fHistMcPid[j] = new TH1F();
131 fHistEsdPid[j] = new TH1F();
132 fHistDataEsdPid[j] = new TH1F();
135 for(Int_t i=0; i<3; i++){
136 fHistMCvertex[i] = new TH1F();
137 fHistESDvertex[i] = new TH1F();
138 fHistDatavertex[i] = new TH1F();
141 for(Int_t iEle=0; iEle<2; iEle++){
142 fHistDataHfePid[iEle] = new TH1F();
147 //________________________________________________________________________
148 AliHFEdca::AliHFEdca(const AliHFEdca &ref):
154 for(Int_t j=0; j<kNParticles; j++){
155 fHistMcPid[j] = ref.fHistMcPid[j];
156 fHistEsdPid[j] = ref.fHistEsdPid[j];
157 fHistDataEsdPid[j] = ref.fHistDataEsdPid[j];
160 for(Int_t i=0; i<3; i++){
161 fHistMCvertex[i] = ref.fHistMCvertex[i];
162 fHistESDvertex[i] = ref.fHistESDvertex[i];
163 fHistDatavertex[i] = ref.fHistDatavertex[i];
166 for(Int_t iEle=0; iEle<2; iEle++){
167 fHistDataHfePid[iEle] = ref.fHistDataHfePid[iEle];
171 //_______________________________________________________________________________________________
172 AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
175 // Assignment operator
178 if(this == &ref) return *this;
179 TObject::operator=(ref);
184 //________________________________________________________________________
185 AliHFEdca::~AliHFEdca()
187 // default destructor
189 for(Int_t j=0; j<kNParticles; j++){
190 for(Int_t i=0; i<kNPtBins; i++){
191 if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
192 if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
194 if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
195 if(fHistDcaZPull[j][i]) delete fHistDcaZPull[j][i];
197 if(fHistDcaXY[j][i]) delete fHistDcaXY[j][i];
198 if(fHistDcaZ[j][i]) delete fHistDcaZ[j][i];
200 if(j<(kNParticles-2)){
201 if(fHistEPDcaXYRes[j][i]) delete fHistEPDcaXYRes[j][i];
202 if(fHistEPDcaZRes[j][i]) delete fHistEPDcaZRes[j][i];
204 if(fHistEPDcaXYPull[j][i]) delete fHistEPDcaXYPull[j][i];
205 if(fHistEPDcaZPull[j][i]) delete fHistEPDcaZPull[j][i];
207 if(fHistEPDcaXY[j][i]) delete fHistEPDcaXY[j][i];
208 if(fHistEPDcaZ[j][i]) delete fHistEPDcaZ[j][i];
211 if(fHistKFDcaXY[j][i]) delete fHistKFDcaXY[j][i];
212 if(fHistKFDcaZ[j][i]) delete fHistKFDcaZ[j][i];
214 if(fHistDataDcaXY[j][i]) delete fHistDataDcaXY[j][i];
215 if(fHistDataDcaZ[j][i]) delete fHistDataDcaZ[j][i];
216 if(fHistDataWoDcaXY[j][i]) delete fHistDataWoDcaXY[j][i];
217 if(fHistDataWoDcaZ[j][i]) delete fHistDataWoDcaZ[j][i];
219 if(fHistDataDcaXYPull[j][i]) delete fHistDataDcaXYPull[j][i];
220 if(fHistDataDcaZPull[j][i]) delete fHistDataDcaZPull[j][i];
221 if(fHistDataWoDcaXYPull[j][i]) delete fHistDataWoDcaXYPull[j][i];
222 if(fHistDataWoDcaZPull[j][i]) delete fHistDataWoDcaZPull[j][i];
225 if(fHistMcPid[j]) delete fHistMcPid[j];
226 if(fHistEsdPid[j]) delete fHistEsdPid[j];
227 if(fHistDataEsdPid[j]) delete fHistDataEsdPid[j];
230 for(Int_t i=0; i<3; i++){
231 if(fHistMCvertex[i]) delete fHistMCvertex[i];
232 if(fHistESDvertex[i]) delete fHistESDvertex[i];
233 if(fHistDatavertex[i]) delete fHistDatavertex[i];
237 for(Int_t iEle=0; iEle<2; iEle++){
238 for(Int_t iPt=0; iPt<kNPtBins; iPt++){
239 if(fHistHPDcaXYRes[iEle][iPt]) delete fHistHPDcaXYRes[iEle][iPt];
240 if(fHistHPDcaZRes[iEle][iPt]) delete fHistHPDcaZRes[iEle][iPt];
241 if(fHistHPDcaXYPull[iEle][iPt]) delete fHistHPDcaXYPull[iEle][iPt];
242 if(fHistHPDcaZPull[iEle][iPt]) delete fHistHPDcaZPull[iEle][iPt];
243 if(fHistHPDcaXY[iEle][iPt]) delete fHistHPDcaXY[iEle][iPt];
244 if(fHistHPDcaZ[iEle][iPt]) delete fHistHPDcaZ[iEle][iPt];
248 if(fHistHPDataDcaXY[iEle][iPt]) delete fHistHPDataDcaXY[iEle][iPt];
249 if(fHistHPDataDcaZ[iEle][iPt]) delete fHistHPDataDcaZ[iEle][iPt];
250 if(fHistHPDataDcaXYPull[iEle][iPt]) delete fHistHPDataDcaXYPull[iEle][iPt];
251 if(fHistHPDataDcaZPull[iEle][iPt]) delete fHistHPDataDcaZPull[iEle][iPt];
254 for(Int_t i=0; i<2; i++)
255 if(fHistHfePid[iEle][i]) delete fHistHfePid[iEle][i];
257 if(fHistDataHfePid[iEle]) delete fHistDataHfePid[iEle];
261 if(fStat) delete fStat;
263 //Printf("analysis done\n");
267 //________________________________________________________________________
268 void AliHFEdca::InitAnalysis(){
270 //Printf("initialize analysis\n");
275 //________________________________________________________________________
276 void AliHFEdca::PostAnalysis() const
279 // moved to dcaPostAnalysis.C
284 //________________________________________________________________________
285 void AliHFEdca::CreateHistogramsResidual(TList *residualList){
290 // fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
291 // fHistDcaZRes[kNParticles][kNPtBins]=0x0;
292 // fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
293 // fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
295 const Int_t nBins = 1000;
296 const Float_t maxXYBin = 1000.;
297 const Float_t maxZBin = 1000.;
300 for(Int_t k=0; k<kNDcaVar; k++){
301 TString histTitle((const char*)fgkResDcaVarTitle[k]);
303 for(Int_t j=0; j<kNParticles; j++){
304 for(Int_t i=0; i<kNPtBins; i++){
306 TString histName((const char*)fgkParticles[j]);
307 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
308 TString histEPName((const char*)fgkParticles[j]);
309 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
312 fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
313 fHistDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);
314 if(j<(kNParticles-2)){
315 fHistEPDcaXYRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
316 fHistEPDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);}
319 fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
320 fHistDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]);
321 if(j<(kNParticles-2)){
322 fHistEPDcaZRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
323 fHistEPDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]); }
329 // TList *fResidualList = 0;
330 residualList->SetOwner();
331 residualList->SetName("residual");
332 for(Int_t iPart=0; iPart<kNParticles; iPart++){
333 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
334 residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
335 residualList->Add(fHistDcaZRes[iPart][iPtBin]);
336 if(iPart<(kNParticles-2)){
337 residualList->Add(fHistEPDcaXYRes[iPart][iPtBin]);
338 residualList->Add(fHistEPDcaZRes[iPart][iPtBin]);
340 } // loop over pt bins
341 } // loop over particles (pos, neg)
348 //________________________________________________________________________
349 void AliHFEdca::CreateHistogramsPull(TList *pullList){
353 const Int_t nBins = 1000;
354 const Float_t maxXYBin = 20.;
355 const Float_t maxZBin = 20.;
358 // for pull -----------------------------------------------------------------------
359 // fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
360 // fHistDcaZPull[kNParticles][kNPtBins]=0x0;
361 // fHistEPDcaXYPull[kNParticles-2][kNPtBins]=0x0;
362 // fHistEPDcaZPull[kNParticles-2][kNPtBins]=0x0;
365 for(Int_t k=0; k<kNDcaVar; k++){
366 TString histTitle((const char*)fgkPullDcaVarTitle[k]);
368 for(Int_t j=0; j<kNParticles; j++){
369 for(Int_t i=0; i<kNPtBins; i++){
371 TString histName((const char*)fgkParticles[j]);
372 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
373 TString histEPName((const char*)fgkParticles[j]);
374 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
377 fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
378 fHistDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
379 if(j<(kNParticles-2)) {
380 fHistEPDcaXYPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
381 fHistEPDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
384 fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
385 fHistDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
386 if(j<(kNParticles-2)) {
387 fHistEPDcaZPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
388 fHistEPDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
394 // TList *fPullList = 0;
395 pullList->SetOwner();
396 pullList->SetName("pull");
397 for(Int_t iPart=0; iPart<kNParticles; iPart++){
398 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
399 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
400 pullList->Add(fHistDcaZPull[iPart][iPtBin]);
401 if(iPart<(kNParticles-2)){
402 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
403 pullList->Add(fHistDcaZPull[iPart][iPtBin]); }
404 } // loop over pt bins
405 } // loop over particles (pos, neg)
410 //________________________________________________________________________
411 void AliHFEdca::CreateHistogramsDca(TList *dcaList){
413 // define histograms: MC dca
418 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
419 fStat->SetMarkerStyle(20);
420 fStat->SetMarkerColor(3);
421 fStat->SetMarkerSize(1);
424 // fHistDcaXY[kNParticles][kNPtBins]=0x0;
425 // fHistDcaZ[kNParticles][kNPtBins]=0x0;
426 // fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
427 // fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
429 const Int_t nBins = 1000;
430 const Float_t maxXYBin = 1000.;
431 const Float_t maxZBin = 1000.;
434 for(Int_t k=0; k<kNDcaVar; k++){
435 TString histTitle((const char*)fgkDcaVarTitle[k]);
437 for(Int_t j=0; j<kNParticles; j++){
438 for(Int_t i=0; i<kNPtBins; i++){
440 TString histName((const char*)fgkParticles[j]);
441 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
443 TString histNameEP((const char*)fgkParticles[j]);
444 histNameEP += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
447 fHistDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
448 fHistDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
450 if(j<(kNParticles-2)){
451 fHistEPDcaXY[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
452 fHistEPDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);}
455 fHistDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
456 fHistDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
457 if(j<(kNParticles-2)){
458 fHistEPDcaZ[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxZBin, maxZBin);
459 fHistEPDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);}
465 // TList *fDcaList = 0;
467 dcaList->SetName("mcDcaDistr");
468 for(Int_t iPart=0; iPart<kNParticles; iPart++){
469 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
470 dcaList->Add(fHistDcaXY[iPart][iPtBin]);
471 dcaList->Add(fHistDcaZ[iPart][iPtBin]);
472 if(iPart<(kNParticles-2)) {
473 dcaList->Add(fHistEPDcaXY[iPart][iPtBin]);
474 dcaList->Add(fHistEPDcaZ[iPart][iPtBin]); }
475 } // loop over pt bins
476 } // loop over particles (pos, neg)
482 //________________________________________________________________________
483 void AliHFEdca::CreateHistogramsKfDca(TList *kfDcaList){
485 // define histograms: MC dca
490 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
491 fStat->SetMarkerStyle(20);
492 fStat->SetMarkerColor(3);
493 fStat->SetMarkerSize(1);
496 // fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
497 // fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
499 const Int_t nBins = 1000;
500 const Float_t maxXYBin = 1000.;
501 const Float_t maxZBin = 1000.;
504 for(Int_t k=0; k<kNDcaVar; k++){
505 TString histTitle((const char*)fgkDcaVarTitle[k]);
507 for(Int_t j=0; j<kNParticles; j++){
508 for(Int_t i=0; i<kNPtBins; i++){
509 TString histNameKF((const char*)fgkParticles[j]);
510 histNameKF += Form("_MCpid_KF%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
513 fHistKFDcaXY[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
514 fHistKFDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
517 fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
518 fHistKFDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
524 kfDcaList->SetOwner();
525 kfDcaList->SetName("mcKfDcaDistr");
526 for(Int_t iPart=0; iPart<kNParticles; iPart++){
527 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
528 kfDcaList->Add(fHistKFDcaXY[iPart][iPtBin]);
529 kfDcaList->Add(fHistKFDcaZ[iPart][iPtBin]);
530 } // loop over pt bins
531 } // loop over particles (pos, neg)
533 kfDcaList->Add(fStat);
538 //________________________________________________________________________
539 void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
541 // define histograms: real Data
545 // fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
546 // fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
547 // fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
548 // fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
550 const Int_t nBins = 1000;
551 const Float_t maxXYBin = 1000.;
552 const Float_t maxZBin = 1000.;
554 for(Int_t k=0; k<kNDcaVar; k++){
555 TString histTitle((const char*)fgkDcaVarTitle[k]);
556 for(Int_t j=0; j<kNParticles; j++){
557 for(Int_t i=0; i<kNPtBins; i++){
559 TString histName((const char*)fgkParticles[j]);
560 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
562 TString histNameWo((const char*)fgkParticles[j]);
563 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
566 fHistDataDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
567 fHistDataDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
569 fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
570 fHistDataWoDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
573 fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
574 fHistDataDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
576 fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
577 fHistDataWoDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
583 // TList *fDcaList = 0;
584 dataDcaList->SetOwner();
585 dataDcaList->SetName("dataDcaDistr");
586 for(Int_t iPart=0; iPart<kNParticles; iPart++){
587 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
588 dataDcaList->Add(fHistDataDcaXY[iPart][iPtBin]);
589 dataDcaList->Add(fHistDataDcaZ[iPart][iPtBin]);
591 dataDcaList->Add(fHistDataWoDcaXY[iPart][iPtBin]);
592 dataDcaList->Add(fHistDataWoDcaZ[iPart][iPtBin]);
593 } // loop over pt bins
594 } // loop over particles (pos, neg)
600 //________________________________________________________________________
601 void AliHFEdca::CreateHistogramsDataPull(TList *dataPullList){
605 const Int_t nBins = 1000;
606 const Float_t maxXYBin = 20.;
607 const Float_t maxZBin = 20.;
609 // for pull -----------------------------------------------------------------------
610 // fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
611 // fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
613 // fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
614 // fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
617 for(Int_t k=0; k<kNDcaVar; k++){
618 TString histTitle((const char*)fgkPullDataDcaVarTitle[k]);
620 for(Int_t j=0; j<kNParticles; j++){
621 for(Int_t i=0; i<kNPtBins; i++){
623 TString histName((const char*)fgkParticles[j]);
624 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
626 TString histNameWo((const char*)fgkParticles[j]);
627 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
630 fHistDataDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
631 fHistDataDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
633 fHistDataWoDcaXYPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
634 fHistDataWoDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
637 fHistDataDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
638 fHistDataDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
640 fHistDataWoDcaZPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
641 fHistDataWoDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
647 // TList *fDataPullList = 0;
648 dataPullList->SetOwner();
649 dataPullList->SetName("dataPull");
650 for(Int_t iPart=0; iPart<kNParticles; iPart++){
651 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
652 dataPullList->Add(fHistDataDcaXYPull[iPart][iPtBin]);
653 dataPullList->Add(fHistDataDcaZPull[iPart][iPtBin]);
655 dataPullList->Add(fHistDataWoDcaXYPull[iPart][iPtBin]);
656 dataPullList->Add(fHistDataWoDcaZPull[iPart][iPtBin]);
657 } // loop over pt bins
658 } // loop over particles (pos, neg)
662 //________________________________________________________________________
663 void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
665 // define histograms: vertex
669 // fHistMCvertex[kNVertexVar]=0x0;
670 // fHistESDvertex[kNVertexVar]=0x0;
672 const Int_t nBins = 1000;
673 const Float_t minXBin = -0.2e4;
674 const Float_t maxXBin = 0.2e4;
675 const Float_t minYBin = -0.5e4;
676 const Float_t maxYBin = 0.5e4;
677 const Float_t minZBin = -1.5e5;
678 const Float_t maxZBin = 1.5e5;
680 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
681 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
683 for(Int_t k=0; k<kNVertexVar; k++){
684 TString histTitle((const char*)fgkVertexVarTitle[k]);
685 TString histNameMC((const char*)fgkVertexVar[k]);
686 histNameMC += Form("_MC");
687 TString histNameESD((const char*)fgkVertexVar[k]);
688 histNameESD += Form("_ESD");
690 fHistMCvertex[k] = new TH1F((const char*)histNameMC, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
691 fHistMCvertex[k]->SetLineColor(k+2);
693 fHistESDvertex[k] = new TH1F((const char*)histNameESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
694 fHistESDvertex[k]->SetLineColor(k+2);
697 vertexList->SetOwner();
698 vertexList->SetName("vertexDistr");
700 for(Int_t k=0; k<kNVertexVar; k++){
701 vertexList->Add(fHistMCvertex[k]);
702 vertexList->Add(fHistESDvertex[k]);
709 //________________________________________________________________________
710 void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
712 // define histograms: vertex
716 // fHistDatavertex[kNVertexVar]=0x0;
718 const Int_t nBins = 1000;
719 const Float_t minXBin = -0.2e4;
720 const Float_t maxXBin = 0.2e4;
721 const Float_t minYBin = -0.5e4;
722 const Float_t maxYBin = 0.5e4;
723 const Float_t minZBin = -1.5e5;
724 const Float_t maxZBin = 1.5e5;
726 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
727 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
729 for(Int_t k=0; k<kNVertexVar; k++){
730 TString histTitle((const char*)fgkVertexVarTitle[k]);
731 TString histNameDataESD((const char*)fgkVertexVar[k]);
732 histNameDataESD += Form("_data");
734 fHistDatavertex[k] = new TH1F((const char*)histNameDataESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
735 fHistDatavertex[k]->SetLineColor(k+2);
738 // TList *fVDaraVertexList = 0;
739 dataVertexList->SetOwner();
740 dataVertexList->SetName("dataVertexDistr");
742 for(Int_t k=0; k<kNVertexVar; k++){
743 dataVertexList->Add(fHistDatavertex[k]);
748 //_______________________________________________________________________________________________
749 void AliHFEdca::CreateHistogramsPid(TList *mcPidList){
751 // define histograms which fills combined PID
754 const Char_t *mcOResd[2]={"mcPt", "esdPt"};
756 for(Int_t iPart=0; iPart<kNParticles; iPart++){
757 TString histTitleMc((const char*)fgkParticles[iPart]);
758 TString histTitleEsd((const char*)fgkParticles[iPart]);
759 histTitleMc += Form("_McPid_%s;p_{T} [GeV/c];counts", mcOResd[0]);
760 histTitleEsd += Form("_EsdPid_%s;p_{T} [GeV/c];counts", mcOResd[1]);
762 TString histNameMc((const char*)fgkParticles[iPart]);
763 TString histNameEsd((const char*)fgkParticles[iPart]);
764 histNameMc+=Form("_McPid_%s", mcOResd[0]);
765 histNameEsd+=Form("_EsdPid_%s", mcOResd[1]);
767 fHistMcPid[iPart] = new TH1F(histNameMc, histTitleMc, kNPtBins, fgkPtIntv);
768 fHistEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
772 mcPidList->SetOwner();
773 mcPidList->SetName("combinedPid");
775 for(Int_t iPart=0; iPart<kNParticles; iPart++){
776 mcPidList->Add(fHistMcPid[iPart]);
777 mcPidList->Add(fHistEsdPid[iPart]);
782 //_______________________________________________________________________________________________
783 void AliHFEdca::CreateHistogramsDataPid(TList *pidList){
785 // define histograms which fills combined PID: data
790 for(Int_t iPart=0; iPart<kNParticles; iPart++){
791 TString histTitleEsd((const char*)fgkParticles[iPart]);
792 histTitleEsd+=Form("_DataEsdPid_esdPt;p_{T} [GeV/c];counts");
793 TString histNameEsd((const char*)fgkParticles[iPart]);
794 histNameEsd+=Form("_DataEsdPid");
796 fHistDataEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
801 pidList->SetName("dataCombinedPid");
803 for(Int_t iPart=0; iPart<kNParticles; iPart++)
804 pidList->Add(fHistDataEsdPid[iPart]);
811 //_______________________________________________________________________________________________
812 void AliHFEdca::FillHistogramsDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
816 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
818 mcPrimV[0] = mcPrimVtx->GetX();
819 mcPrimV[1] = mcPrimVtx->GetY();
820 mcPrimV[2] = mcPrimVtx->GetZ();
822 Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
824 // filling historgams track by track
825 // obtaining reconstructed dca ------------------------------------------------------------------
827 Float_t esdpx = track->Px();
828 Float_t esdpy = track->Py();
829 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
831 // obtaining errors of dca ------------------------------------------------------------------
832 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
834 primV[0] = primVtx->GetXv();
835 primV[1] = primVtx->GetYv();
836 primV[2] = primVtx->GetZv();
838 Float_t magneticField = 0; // initialized as 5kG
839 magneticField = esdEvent->GetMagneticField(); // in kG
841 Double_t beampiperadius=3.;
842 Double_t dz[2]; // error of dca in cm
845 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
847 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
849 TParticle *part = mctrack->Particle();
851 Float_t vx = part->Vx(); // in cm
852 Float_t vy = part->Vy(); // in cm
853 Float_t vz = part->Vz(); // in cm
855 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
857 Float_t mcpx = part->Px();
858 Float_t mcpy = part->Py();
859 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
861 Int_t pdg = part->GetPdgCode();
862 Int_t esdPid = GetCombinedPid(track);
865 if(pdg==kPDGelectron || pdg==kPDGmuon
866 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
868 // calculate mcDca ------------------------------------------------------------------
869 const Float_t conv[2] = {1.783/1.6, 2.99792458};
870 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
872 Float_t nx = esdpx/mcpt;
873 Float_t ny = esdpy/mcpt;
876 radius = TMath::Abs(radiusMc);
878 Double_t dxy = vxy - mcVtxXY; // in cm
879 Double_t dvx = vx - mcPrimV[0]; // in cm
880 Double_t dvy = vy - mcPrimV[1]; // in cm
882 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
884 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
885 Double_t residual[2] = {0, 0};
886 Double_t pull[2] = {0, 0};
887 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
888 for(Int_t i=0; i<2; i++){
889 residual[i] = dz[i] - mcDca[i]; // in centimeters
890 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
894 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
896 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
897 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
898 if(pdg==fgkPdgParticle[iPart]) {
899 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
900 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
901 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
902 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
903 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
904 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
907 if(esdPid==fgkPdgParticle[iPart]) {
908 fHistEPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
909 fHistEPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
910 fHistEPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
911 fHistEPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
912 fHistEPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
913 fHistEPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
921 } // particle id loop
923 // for charged particles: no pid
924 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
925 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
927 if(charge>0) iPart = 11;
928 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
929 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
930 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
931 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
932 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
933 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
941 //_______________________________________________________________________________________________
942 void AliHFEdca::FillHistogramsKfDca(AliESDEvent * const esdEvent, AliESDtrack * const track, const AliMCEvent * const mcEvent)
946 // filling historgams track by track
948 // obtaining reconstructed dca ------------------------------------------------------------------
949 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
950 Float_t magneticField = 0; // initialized as 5kG
951 magneticField = esdEvent->GetMagneticField(); // in kG
953 Float_t esdpx = track->Px();
954 Float_t esdpy = track->Py();
955 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
957 Int_t charge = track->Charge();
959 Double_t beampiperadius=3.;
960 Double_t dz[2]; // error of dca in cm
962 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
964 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
966 TParticle *part = mctrack->Particle();
967 Int_t pdg = part->GetPdgCode();
969 // calculate dca using AliKFParticle class------------------------------------------------------------------
970 Double_t kfdz[3] = {0, 0, 0};
971 Double_t kfdzwith[3] = {0, 0, 0};
973 Int_t trkID = track->GetID();
975 AliKFParticle::SetField(magneticField);
976 AliKFParticle kfParticle(*track, pdg);
978 // prepare kfprimary vertex
979 AliKFVertex kfESDprimary;
980 // Reconstruct Primary Vertex (with ESD tracks)
981 Int_t n=primVtx->GetNIndices();
982 if (n>0 && primVtx->GetStatus()){
983 kfESDprimary = AliKFVertex(*primVtx);
985 Double_t dcaXYWithTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
986 Double_t dcaWithTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
987 Double_t dcaZWithTrk = 0;
988 if(TMath::Abs(dcaWithTrk)>=TMath::Abs(dcaXYWithTrk))
989 dcaZWithTrk =TMath::Sqrt(dcaWithTrk*dcaWithTrk-dcaXYWithTrk*dcaXYWithTrk)*((dz[1]*-1<=0)?1:-1);
990 kfdzwith[0] = dcaXYWithTrk;
991 kfdzwith[1] = dcaZWithTrk;
992 kfdzwith[2] = dcaWithTrk; // with current track
994 Double_t dcaXYWoTrk = 0;
995 Double_t dcaZWoTrk = 0;
996 Double_t dcaWoTrk = 0;
998 UShort_t *priIndex = primVtx->GetIndices();
1000 for (Int_t i=0;i<n;i++){
1002 Int_t idx = Int_t(priIndex[i]);
1004 kfESDprimary -= kfParticle;
1005 dcaXYWoTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
1006 dcaWoTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
1007 if((dcaWoTrk-dcaXYWoTrk)>=0)
1008 dcaZWoTrk = TMath::Abs(dcaWoTrk*dcaWoTrk - dcaXYWoTrk*dcaXYWoTrk)*((dz[1]*-1<=0)?1:-1);
1009 } // remove current track from this calculation
1010 } // loop over all primary vertex contributors
1013 kfdz[0] = dcaXYWoTrk;
1014 kfdz[1] = dcaZWoTrk;
1017 } // only if n contributor > 0 and primVtx constructed
1021 if(dz[0]!=0 && dz[0]*kfdzwith[0]>0 && TMath::Abs(kfdzwith[0]/dz[0])>0.9999 && TMath::Abs(kfdzwith[0]/dz[0])<1.0001)fStat->Fill(1);; // same
1022 if(dz[0]!=0 && dz[0]*kfdzwith[0]<0 && TMath::Abs(kfdzwith[0]/dz[0])>0.9999 && TMath::Abs(kfdzwith[0]/dz[0])<1.0001) fStat->Fill(2); // swapped sign
1023 if(kfdzwith[0]==0 && dz[0]!=0) fStat->Fill(3); // 0 from KF particle (with current track)
1025 if(dz[0]!=0 && dz[0]*kfdz[0]>0 && TMath::Abs(kfdz[0]/dz[0])>0.8 && TMath::Abs(kfdz[0]/dz[0])<1.2) fStat->Fill(-1);; // same
1026 if(dz[0]!=0 && dz[0]*kfdz[0]<0 && TMath::Abs(kfdz[0]/dz[0])>0.8 && TMath::Abs(kfdz[0]/dz[0]) <1.2) fStat->Fill(-2); // swapped sign
1027 if(kfdz[0]==0 && dz[0]!=0) fStat->Fill(-3); // 0 from KF particle (without current track)
1029 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1031 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1032 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1033 if(pdg==fgkPdgParticle[iPart]) {
1034 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1.0e4);
1035 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1.0e4);
1042 } // particle id loop
1044 // for charged particles: no pid
1045 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1046 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1048 if(charge>0) iPart = 11;
1049 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1050 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1060 //_______________________________________________________________________________________________
1061 void AliHFEdca::FillHistogramsVtx(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
1065 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1066 Double_t mcPrimV[3];
1067 mcPrimV[0] = mcPrimVtx->GetX();
1068 mcPrimV[1] = mcPrimVtx->GetY();
1069 mcPrimV[2] = mcPrimVtx->GetZ();
1071 // obtaining errors of dca ------------------------------------------------------------------
1072 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1074 primV[0] = primVtx->GetXv();
1075 primV[1] = primVtx->GetYv();
1076 primV[2] = primVtx->GetZv();
1078 for(Int_t i=0; i<kNVertexVar; i++){
1079 fHistMCvertex[i]->Fill(mcPrimV[i]*1.0e4);
1080 fHistESDvertex[i]->Fill(primV[i]*1.0e4);
1085 //_______________________________________________________________________________________________
1086 void AliHFEdca::FillHistogramsPid(AliESDtrack * const track, const AliMCEvent * const mcEvent)
1090 // filling historgams track by track
1091 Float_t esdpx = track->Px();
1092 Float_t esdpy = track->Py();
1093 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1095 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1096 if(!mctrack) return;
1097 TParticle *part = mctrack->Particle();
1099 Float_t mcpx = part->Px();
1100 Float_t mcpy = part->Py();
1101 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1103 Int_t pdg = part->GetPdgCode();
1104 Int_t esdPid = GetCombinedPid(track);
1107 Double_t ptMom[2] = {mcpt, esdpt};
1109 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1110 if(pdg==fgkPdgParticle[iPart]) // pid all by MC
1111 fHistMcPid[iPart]->Fill(ptMom[0]);
1113 if(esdPid==fgkPdgParticle[iPart]) // pid all by combined pid
1114 fHistEsdPid[iPart]->Fill(ptMom[1]);
1115 } // loop over particles
1118 if(pdg==kPDGelectron || pdg==kPDGmuon || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton)
1119 fHistMcPid[10]->Fill(ptMom[0]);
1120 if(pdg==-kPDGelectron || pdg==-kPDGmuon || pdg==kPDGpion || pdg==kPDGkaon || pdg==kPDGproton)
1121 fHistMcPid[11]->Fill(ptMom[0]);
1122 if(esdPid==kPDGelectron || esdPid==kPDGmuon || esdPid==-kPDGpion || esdPid==-kPDGkaon || esdPid==-kPDGproton)
1123 fHistEsdPid[10]->Fill(ptMom[1]);
1124 if(esdPid==-kPDGelectron || esdPid==-kPDGmuon || esdPid==kPDGpion || esdPid==kPDGkaon || esdPid==kPDGproton)
1125 fHistEsdPid[11]->Fill(ptMom[1]);
1129 ////_______________________________________________________________________________________________
1130 void AliHFEdca::FillHistogramsDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1132 // filling historgams track by track
1133 // obtaining reconstructed dca --------------------------------------------------------------
1135 Float_t esdpx = track->Px();
1136 Float_t esdpy = track->Py();
1137 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1138 Int_t charge = track->Charge();
1140 // obtaining errors of dca ------------------------------------------------------------------
1141 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1143 primV[0] = primVtx->GetXv();
1144 primV[1] = primVtx->GetYv();
1145 primV[2] = primVtx->GetZv();
1148 Float_t magneticField = 0; // initialized as 5kG
1149 magneticField = esdEvent->GetMagneticField(); // in kG
1151 Double_t beampiperadius=3.;
1152 Double_t dz[2]; // error of dca in cm
1153 Double_t covardz[3];
1155 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1158 Double_t pull[2] = {0, 0};
1159 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1160 for(Int_t i=0; i<2; i++){
1161 if(error[i]!=0)pull[i] = dz[i]/error[i]; // unitless
1164 // get dca when current track is not included
1166 Double_t dzwo[2], covardzwo[3];
1167 Double_t pullwo[2] = {0, 0};
1168 if(!track->PropagateToDCA(vtxESDSkip, magneticField, beampiperadius, dzwo, covardzwo)) return; // protection
1170 Double_t errorwo[2] ={TMath::Sqrt(TMath::Abs(covardzwo[0])), TMath::Sqrt(TMath::Abs(covardzwo[2]))};
1171 for(Int_t i=0; i<2; i++){
1172 if(errorwo[i]!=0) pullwo[i] = dzwo[i]/errorwo[i]; // unitless
1176 Int_t esdPid = GetCombinedPid(track);
1178 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1180 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1181 if(esdPid==fgkPdgParticle[iPart] && (esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1])) {
1182 fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1183 fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1184 fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1185 fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1186 // w/o current track
1187 fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1188 fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1189 fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1190 fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1197 // for charged particles
1198 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1199 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1201 if(charge>0) iPart = 11;
1202 fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1203 fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1204 fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1205 fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1206 // without current track
1207 fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1208 fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1209 fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1210 fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1219 //_______________________________________________________________________________________________
1220 void AliHFEdca::FillHistogramsDataVtx(AliESDEvent * const esdEvent)
1224 // obtaining errors of dca ------------------------------------------------------------------
1225 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1227 primV[0] = primVtx->GetXv();
1228 primV[1] = primVtx->GetYv();
1229 primV[2] = primVtx->GetZv();
1231 // require events with at least 3 contributors for primary vertex construction
1232 Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1233 if(nEsdPrimVtxCtrb<1) return; // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1234 for(Int_t i=0; i<kNVertexVar; i++)
1235 fHistDatavertex[i]->Fill(primV[i]*1e4);
1240 ////_______________________________________________________________________________________________
1241 void AliHFEdca::FillHistogramsDataPid(AliESDtrack * const track)
1243 // filling historgams track by track
1244 // obtaining reconstructed dca --------------------------------------------------------------
1246 Float_t esdpx = track->Px();
1247 Float_t esdpy = track->Py();
1248 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1249 Int_t charge = track->Charge();
1251 Int_t esdPid = GetCombinedPid(track);
1253 for(Int_t iPart=0; iPart<kNParticles; iPart++){
1254 if(iPart<(kNParticles-2)){
1255 if(esdPid==fgkPdgParticle[iPart]) fHistDataEsdPid[iPart]->Fill(esdpt);
1259 if(charge<0) fHistDataEsdPid[10]->Fill(esdpt);
1260 if(charge>0) fHistDataEsdPid[11]->Fill(esdpt);
1265 //_________________________________________________________________________________________________
1266 void AliHFEdca::ApplyExtraCuts(AliESDEvent * const esdEvent, Int_t nMinPrimVtxContributor)
1270 // only one extra cut, number of contributors to each primary vertex
1273 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1274 Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1275 if(nEsdPrimVtxCtrb<nMinPrimVtxContributor) return;
1276 // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1280 //_____________________________________________________
1281 Int_t AliHFEdca::GetCombinedPid(AliESDtrack *const track)
1284 // combined detector pid
1285 Double_t prob[AliPID::kSPECIES];
1286 track->GetESDpid(prob);
1289 Double_t priors[AliPID::kSPECIESN];
1296 Int_t charge = track->Charge();
1300 pid.SetPriors(priors);
1301 pid.SetProbabilities(prob);
1303 // identify particle as the most probable
1305 Double_t pelectron = pid.GetProbability(AliPID::kElectron);
1306 if(pelectron > pid.GetProbability(AliPID::kMuon) &&
1307 pelectron > pid.GetProbability(AliPID::kPion) &&
1308 pelectron > pid.GetProbability(AliPID::kKaon) &&
1309 pelectron > pid.GetProbability(AliPID::kProton) ) esdPid = -kPDGelectron;
1311 Double_t pmuon = pid.GetProbability(AliPID::kMuon);
1312 if(pmuon > pid.GetProbability(AliPID::kElectron) &&
1313 pmuon > pid.GetProbability(AliPID::kPion) &&
1314 pmuon > pid.GetProbability(AliPID::kKaon) &&
1315 pmuon > pid.GetProbability(AliPID::kProton) ) esdPid = -kPDGmuon;
1317 Double_t ppion = pid.GetProbability(AliPID::kPion);
1318 if(ppion > pid.GetProbability(AliPID::kElectron) &&
1319 ppion > pid.GetProbability(AliPID::kMuon) &&
1320 ppion > pid.GetProbability(AliPID::kKaon) &&
1321 ppion > pid.GetProbability(AliPID::kProton) ) esdPid = kPDGpion;
1323 Double_t pkaon = pid.GetProbability(AliPID::kKaon);
1324 if(pkaon > pid.GetProbability(AliPID::kElectron) &&
1325 pkaon > pid.GetProbability(AliPID::kMuon) &&
1326 pkaon > pid.GetProbability(AliPID::kPion) &&
1327 pkaon > pid.GetProbability(AliPID::kProton) ) esdPid = kPDGkaon;
1329 Double_t pproton = pid.GetProbability(AliPID::kProton);
1330 if(pproton > pid.GetProbability(AliPID::kElectron) &&
1331 pproton > pid.GetProbability(AliPID::kMuon) &&
1332 pproton > pid.GetProbability(AliPID::kPion) &&
1333 pproton > pid.GetProbability(AliPID::kKaon) ) esdPid = kPDGproton;
1336 return charge*esdPid;
1343 //________________________________________________________________________
1344 void AliHFEdca::CreateHistogramsHfeDca(TList *hfeDcaList){
1346 // define histograms: hfe pid electrons in MC
1349 const Int_t nBinsDca = 1000;
1350 const Float_t maxXYBinDca = 1000.;
1351 const Float_t maxZBinDca = 1000.;
1353 const Int_t nBinsPull = 1000;
1354 const Float_t maxXYBinPull = 20.;
1355 const Float_t maxZBinPull = 20.;
1357 const Char_t *mcOResd[2]={"mcPt", "esdPt"};
1359 // fHistHPDcaXY[2][kNPtBins]=0x0;
1360 // fHistHPDcaZ[2][kNPtBins]=0x0;
1361 // fHistHPDcaXYRes[2][kNPtBins]=0x0;
1362 // fHistHPDcaZRes[2][kNPtBins]=0x0;
1363 // fHistHPDcaXYPull[2][kNPtBins]=0x0;
1364 // fHistHPDcaZPull[2][kNPtBins]=0x0;
1367 for(Int_t k=0; k<kNDcaVar; k++){
1369 TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1370 TString histTitleRes((const char*)fgkPullDcaVarTitle[k]);
1371 TString histTitlePull((const char*)fgkResDcaVarTitle[k]);
1373 for(Int_t iPart=0; iPart<2; iPart++){
1374 for(Int_t i=0; i<kNPtBins; i++){
1375 TString histHPName((const char*)fgkParticles[iPart*5]);
1376 histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1377 TString histHPNameRes((const char*)fgkParticles[iPart*5]);
1378 histHPNameRes += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1379 TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1380 histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1383 fHistHPDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1384 fHistHPDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1385 fHistHPDcaXYRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1386 fHistHPDcaXYRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1387 fHistHPDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1388 fHistHPDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1392 fHistHPDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1393 fHistHPDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1394 fHistHPDcaZRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1395 fHistHPDcaZRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1396 fHistHPDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxZBinPull, 1+maxZBinPull);
1397 fHistHPDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1403 // fHistHfePid[2][2] = 0x0; //! HFE pid
1404 for(Int_t id=0; id<2; id++){
1405 for(Int_t iPart=0; iPart<2; iPart++){
1406 TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1407 histTitleHfe+=Form("_MC_HfePid_esdPt;p_{T} [GeV/c];counts");
1408 TString histNameHfe((const char*)fgkParticles[iPart*5]);
1409 histNameHfe+=Form("_MC_HfePid_%s", mcOResd[id]);
1410 fHistHfePid[id][iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1414 hfeDcaList->SetOwner();
1415 hfeDcaList->SetName("hfeDca");
1416 for(Int_t iPart=0; iPart<2; iPart++){
1417 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1418 hfeDcaList->Add(fHistHPDcaXY[iPart][iPtBin]);
1419 hfeDcaList->Add(fHistHPDcaZ[iPart][iPtBin]);
1420 hfeDcaList->Add(fHistHPDcaXYRes[iPart][iPtBin]);
1421 hfeDcaList->Add(fHistHPDcaZRes[iPart][iPtBin]);
1422 hfeDcaList->Add(fHistHPDcaXYPull[iPart][iPtBin]);
1423 hfeDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1425 for(Int_t id=0; id<2; id++)
1426 hfeDcaList->Add(fHistHfePid[id][iPart]);
1432 //________________________________________________________________________
1433 void AliHFEdca::CreateHistogramsHfeDataDca(TList *hfeDataDcaList){
1435 // define histograms: hfe pid electrons in data
1438 const Int_t nBinsDca = 1000;
1439 const Float_t maxXYBinDca = 1000.;
1440 const Float_t maxZBinDca = 1000.;
1442 const Int_t nBinsPull = 1000;
1443 const Float_t maxXYBinPull = 20.;
1444 const Float_t maxZBinPull = 20.;
1447 // fHistHPDataDcaXY[2][kNPtBins]=0x0;
1448 // fHistHPDataDcaZ[2][kNPtBins]=0x0;
1449 // fHistHPDataDcaXYPull[2][kNPtBins]=0x0;
1450 // fHistHPDataDcaZPull[2][kNPtBins]=0x0;
1452 for(Int_t k=0; k<kNDcaVar; k++){
1453 TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1454 TString histTitlePull((const char*)fgkPullDcaVarTitle[k]);
1455 for(Int_t iPart=0; iPart<2; iPart++){
1456 for(Int_t i=0; i<kNPtBins; i++){
1457 TString histHPName((const char*)fgkParticles[iPart*5]);
1458 histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1459 TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1460 histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1463 fHistHPDataDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1464 fHistHPDataDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1465 fHistHPDataDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1466 fHistHPDataDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1470 fHistHPDataDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1471 fHistHPDataDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1472 fHistHPDataDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsDca, 1-maxZBinPull, 1+maxZBinPull);
1473 fHistHPDataDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1477 } // 2 particle type
1480 //fHistDataHfePid[2] = 0x0; //! HFE pid
1481 for(Int_t iPart=0; iPart<2; iPart++){
1482 TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1483 histTitleHfe+=Form("_Data_HfePid_esdPt;p_{T} [GeV/c];counts");
1484 TString histNameHfe((const char*)fgkParticles[iPart*5]);
1485 histNameHfe+=Form("_Data_HfePid");
1486 fHistDataHfePid[iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1490 hfeDataDcaList->SetOwner();
1491 hfeDataDcaList->SetName("hfeDataDca");
1492 for(Int_t iPart=0; iPart<2; iPart++){
1493 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1494 hfeDataDcaList->Add(fHistHPDataDcaXY[iPart][iPtBin]);
1495 hfeDataDcaList->Add(fHistHPDataDcaZ[iPart][iPtBin]);
1496 hfeDataDcaList->Add(fHistHPDataDcaXYPull[iPart][iPtBin]);
1497 hfeDataDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1499 hfeDataDcaList->Add(fHistDataHfePid[iPart]);
1507 //_______________________________________________________________________________________________
1508 void AliHFEdca::FillHistogramsHfeDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
1510 // the kHFEpid plugin
1512 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1513 Double_t mcPrimV[3];
1514 mcPrimV[0] = mcPrimVtx->GetX();
1515 mcPrimV[1] = mcPrimVtx->GetY();
1516 mcPrimV[2] = mcPrimVtx->GetZ();
1518 Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
1520 // filling historgams track by track
1521 // obtaining reconstructed dca ------------------------------------------------------------------
1522 Float_t esdpx = track->Px();
1523 Float_t esdpy = track->Py();
1524 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1526 // obtaining errors of dca ------------------------------------------------------------------
1527 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1529 primV[0] = primVtx->GetXv();
1530 primV[1] = primVtx->GetYv();
1531 primV[2] = primVtx->GetZv();
1533 Float_t magneticField = 0; // initialized as 5kG
1534 magneticField = esdEvent->GetMagneticField(); // in kG
1536 Double_t beampiperadius=3.;
1537 Double_t dz[2]; // error of dca in cm
1538 Double_t covardz[3];
1539 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1541 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1542 if(!mctrack) return;
1543 TParticle *part = mctrack->Particle();
1545 Float_t vx = part->Vx(); // in cm
1546 Float_t vy = part->Vy(); // in cm
1547 Float_t vz = part->Vz(); // in cm
1549 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1551 Float_t mcpx = part->Px();
1552 Float_t mcpy = part->Py();
1553 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1555 Int_t pdg = part->GetPdgCode();
1558 if(pdg==kPDGelectron || pdg==kPDGmuon
1559 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1561 // calculate mcDca ------------------------------------------------------------------
1562 const Float_t conv[2] = {1.783/1.6, 2.99792458};
1563 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
1565 Float_t nx = esdpx/mcpt;
1566 Float_t ny = esdpy/mcpt;
1569 radius = TMath::Abs(radiusMc);
1571 Double_t dxy = vxy - mcVtxXY; // in cm
1572 Double_t dvx = vx - mcPrimV[0]; // in cm
1573 Double_t dvy = vy - mcPrimV[1]; // in cm
1575 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
1577 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
1578 Double_t residual[2] = {0, 0};
1579 Double_t pull[2] = {0, 0};
1580 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1581 for(Int_t i=0; i<2; i++){
1582 residual[i] = dz[i] - mcDca[i]; // in centimeters
1583 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
1587 if(track->Charge()<0) iPart = 0; // electron
1588 if(track->Charge()>0) iPart = 1; // positron
1589 if(track->Charge()==0) {
1590 printf("this is not an electron! Check HFEpid method");
1593 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1594 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1595 fHistHPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
1596 fHistHPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
1597 fHistHPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1598 fHistHPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1599 fHistHPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
1600 fHistHPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
1608 fHistHfePid[iPart][0]->Fill(esdpt);
1609 fHistHfePid[iPart][1]->Fill(mcpt);
1614 //_______________________________________________________________________________________________
1615 void AliHFEdca::FillHistogramsHfeDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1617 // filling historgams track by track
1618 // obtaining reconstructed dca --------------------------------------------------------------
1620 Float_t esdpx = track->Px();
1621 Float_t esdpy = track->Py();
1622 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1623 Int_t charge = track->Charge();
1625 // obtaining errors of dca ------------------------------------------------------------------
1626 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!
1628 primV[0] = primVtx->GetXv();
1629 primV[1] = primVtx->GetYv();
1630 primV[2] = primVtx->GetZv();
1632 Float_t magneticField = 0; // initialized as 5kG
1633 magneticField = esdEvent->GetMagneticField(); // in kG
1634 Double_t beampiperadius=3.;
1636 Double_t dz[2]; // error of dca in cm
1637 Double_t covardz[3];
1639 if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
1641 Double_t pull[2] = {0, 0};
1642 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1643 for(Int_t i=0; i<2; i++){
1644 if(error[i]!=0) pull[i] = dz[i]/error[i]; // unitless
1648 if(charge<0) iPart = 0; // electron
1649 if(charge>0) iPart = 1; // positron
1651 printf("this is not an electron! Check HFEpid method\n");
1655 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1656 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]) {
1657 fHistHPDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1658 fHistHPDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1659 fHistHPDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1660 fHistHPDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1666 fHistDataHfePid[iPart]->Fill(esdpt);