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 **************************************************************************/
16 // Class for impact parameter (DCA) of charged particles
17 // Study resolution and pull: prepare for beauty study
20 // Hongyan Yang <hongyan@physi.uni-heidelberg.de>
21 // Carlo Bombonati <carlo.bombonati@cern.ch>
27 #include <TParticle.h>
28 #include "AliMCParticle.h"
29 #include "AliESDtrack.h"
30 #include "AliESDEvent.h"
31 #include "AliMCEvent.h"
32 #include "AliMCVertex.h"
34 #include "AliKFParticle.h"
35 #include "AliKFVertex.h"
37 #include "AliESDVertex.h"
41 #include "AliHFEdca.h"
46 //________________________________________________________________________
47 const Char_t* AliHFEdca::fgkParticles[12] = {
49 "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
50 "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
51 "allNegative", "allPositive"
54 const Int_t AliHFEdca::fgkPdgParticle[10] = {
55 // 11, 13, -211, -233, -2122,
56 // -11, -13, 211, 233, 2122};
57 kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
58 -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
60 //________________________________________________________________________
61 const Int_t AliHFEdca::fgkColorPart[12] = {
62 // colors assigned to particles
63 kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
64 kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
68 //________________________________________________________________________
69 const Float_t AliHFEdca::fgkPtIntv[51] = {
71 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55,
72 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50,
73 1.70, 1.90, 2.10, 2.30, 2.50, 2.70, 2.90, 3.10, 3.30, 3.50,
74 3.80, 4.10, 4.40, 4.70, 5.00, 5.30, 5.60, 5.90, 6.20, 6.50,
75 7.00, 7.50, 8.00, 9.00, 10.0, 11.0, 12.0, 13.0, 15.0, 18.0,
78 //________________________________________________________________________
79 const Char_t *AliHFEdca::fgkDcaVar[2] = {
82 //________________________________________________________________________
83 const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
84 ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
86 //________________________________________________________________________
87 const Char_t *AliHFEdca::fgkVertexVar[3] = {
88 "VertexX", "VertexY", "VertexZ"};
90 //________________________________________________________________________
91 const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
92 ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
94 //________________________________________________________________________
95 const Char_t *AliHFEdca::fgkResDcaVar[2] = {
96 "deltaDcaXY", "deltaDcaZ"};
98 //________________________________________________________________________
99 const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
100 ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
102 //________________________________________________________________________
103 const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
104 "pullDcaXY", "pullDcaZ"
107 //________________________________________________________________________
108 const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
109 ";residual dca_{xy}/(error dca_{xy});counts",
110 ";residual dca_{z}/(error dca_{z});counts"
113 //________________________________________________________________________
114 const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
115 ";dca_{xy}^{data}/error dca_{xy};counts",
116 ";dca_{z}^{data}/error dca_{z};counts"
119 //________________________________________________________________________
120 AliHFEdca::AliHFEdca():
124 // default constructor
126 for(Int_t j=0; j<kNParticles; j++){
127 fHistMcPid[j] = new TH1F();
128 fHistEsdPid[j] = new TH1F();
129 fHistDataEsdPid[j] = new TH1F();
132 for(Int_t i=0; i<3; i++){
133 fHistMCvertex[i] = new TH1F();
134 fHistESDvertex[i] = new TH1F();
135 fHistDatavertex[i] = new TH1F();
138 for(Int_t iEle=0; iEle<2; iEle++){
139 fHistDataHfePid[iEle] = new TH1F();
144 //________________________________________________________________________
145 AliHFEdca::AliHFEdca(const AliHFEdca &ref):
151 for(Int_t j=0; j<kNParticles; j++){
152 fHistMcPid[j] = ref.fHistMcPid[j];
153 fHistEsdPid[j] = ref.fHistEsdPid[j];
154 fHistDataEsdPid[j] = ref.fHistDataEsdPid[j];
157 for(Int_t i=0; i<3; i++){
158 fHistMCvertex[i] = ref.fHistMCvertex[i];
159 fHistESDvertex[i] = ref.fHistESDvertex[i];
160 fHistDatavertex[i] = ref.fHistDatavertex[i];
163 for(Int_t iEle=0; iEle<2; iEle++){
164 fHistDataHfePid[iEle] = ref.fHistDataHfePid[iEle];
168 //_______________________________________________________________________________________________
169 AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
172 // Assignment operator
175 if(this == &ref) return *this;
176 TObject::operator=(ref);
181 //________________________________________________________________________
182 AliHFEdca::~AliHFEdca()
184 // default destructor
186 for(Int_t j=0; j<kNParticles; j++){
187 for(Int_t i=0; i<kNPtBins; i++){
188 if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
189 if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
191 if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
192 if(fHistDcaZPull[j][i]) delete fHistDcaZPull[j][i];
194 if(fHistDcaXY[j][i]) delete fHistDcaXY[j][i];
195 if(fHistDcaZ[j][i]) delete fHistDcaZ[j][i];
197 if(j<(kNParticles-2)){
198 if(fHistEPDcaXYRes[j][i]) delete fHistEPDcaXYRes[j][i];
199 if(fHistEPDcaZRes[j][i]) delete fHistEPDcaZRes[j][i];
201 if(fHistEPDcaXYPull[j][i]) delete fHistEPDcaXYPull[j][i];
202 if(fHistEPDcaZPull[j][i]) delete fHistEPDcaZPull[j][i];
204 if(fHistEPDcaXY[j][i]) delete fHistEPDcaXY[j][i];
205 if(fHistEPDcaZ[j][i]) delete fHistEPDcaZ[j][i];
208 if(fHistKFDcaXY[j][i]) delete fHistKFDcaXY[j][i];
209 if(fHistKFDcaZ[j][i]) delete fHistKFDcaZ[j][i];
211 if(fHistDataDcaXY[j][i]) delete fHistDataDcaXY[j][i];
212 if(fHistDataDcaZ[j][i]) delete fHistDataDcaZ[j][i];
213 if(fHistDataWoDcaXY[j][i]) delete fHistDataWoDcaXY[j][i];
214 if(fHistDataWoDcaZ[j][i]) delete fHistDataWoDcaZ[j][i];
216 if(fHistDataDcaXYPull[j][i]) delete fHistDataDcaXYPull[j][i];
217 if(fHistDataDcaZPull[j][i]) delete fHistDataDcaZPull[j][i];
218 if(fHistDataWoDcaXYPull[j][i]) delete fHistDataWoDcaXYPull[j][i];
219 if(fHistDataWoDcaZPull[j][i]) delete fHistDataWoDcaZPull[j][i];
222 if(fHistMcPid[j]) delete fHistMcPid[j];
223 if(fHistEsdPid[j]) delete fHistEsdPid[j];
224 if(fHistDataEsdPid[j]) delete fHistDataEsdPid[j];
227 for(Int_t i=0; i<3; i++){
228 if(fHistMCvertex[i]) delete fHistMCvertex[i];
229 if(fHistESDvertex[i]) delete fHistESDvertex[i];
230 if(fHistDatavertex[i]) delete fHistDatavertex[i];
234 for(Int_t iEle=0; iEle<2; iEle++){
235 for(Int_t iPt=0; iPt<kNPtBins; iPt++){
236 if(fHistHPDcaXYRes[iEle][iPt]) delete fHistHPDcaXYRes[iEle][iPt];
237 if(fHistHPDcaZRes[iEle][iPt]) delete fHistHPDcaZRes[iEle][iPt];
238 if(fHistHPDcaXYPull[iEle][iPt]) delete fHistHPDcaXYPull[iEle][iPt];
239 if(fHistHPDcaZPull[iEle][iPt]) delete fHistHPDcaZPull[iEle][iPt];
240 if(fHistHPDcaXY[iEle][iPt]) delete fHistHPDcaXY[iEle][iPt];
241 if(fHistHPDcaZ[iEle][iPt]) delete fHistHPDcaZ[iEle][iPt];
245 if(fHistHPDataDcaXY[iEle][iPt]) delete fHistHPDataDcaXY[iEle][iPt];
246 if(fHistHPDataDcaZ[iEle][iPt]) delete fHistHPDataDcaZ[iEle][iPt];
247 if(fHistHPDataDcaXYPull[iEle][iPt]) delete fHistHPDataDcaXYPull[iEle][iPt];
248 if(fHistHPDataDcaZPull[iEle][iPt]) delete fHistHPDataDcaZPull[iEle][iPt];
251 for(Int_t i=0; i<2; i++)
252 if(fHistHfePid[iEle][i]) delete fHistHfePid[iEle][i];
254 if(fHistDataHfePid[iEle]) delete fHistDataHfePid[iEle];
258 if(fStat) delete fStat;
260 //Printf("analysis done\n");
264 //________________________________________________________________________
265 void AliHFEdca::InitAnalysis(){
267 //Printf("initialize analysis\n");
272 //________________________________________________________________________
273 void AliHFEdca::PostAnalysis() const
276 // moved to dcaPostAnalysis.C
281 //________________________________________________________________________
282 void AliHFEdca::CreateHistogramsResidual(TList *residualList){
287 // fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
288 // fHistDcaZRes[kNParticles][kNPtBins]=0x0;
289 // fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
290 // fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
292 const Int_t nBins = 1000;
293 const Float_t maxXYBin = 1000.;
294 const Float_t maxZBin = 1000.;
297 for(Int_t k=0; k<kNDcaVar; k++){
298 TString histTitle((const char*)fgkResDcaVarTitle[k]);
300 for(Int_t j=0; j<kNParticles; j++){
301 for(Int_t i=0; i<kNPtBins; i++){
303 TString histName((const char*)fgkParticles[j]);
304 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
305 TString histEPName((const char*)fgkParticles[j]);
306 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
309 fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
310 fHistDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);
311 if(j<(kNParticles-2)){
312 fHistEPDcaXYRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
313 fHistEPDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);}
316 fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
317 fHistDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]);
318 if(j<(kNParticles-2)){
319 fHistEPDcaZRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
320 fHistEPDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]); }
326 // TList *fResidualList = 0;
327 residualList->SetOwner();
328 residualList->SetName("residual");
329 for(Int_t iPart=0; iPart<kNParticles; iPart++){
330 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
331 residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
332 residualList->Add(fHistDcaZRes[iPart][iPtBin]);
333 if(iPart<(kNParticles-2)){
334 residualList->Add(fHistEPDcaXYRes[iPart][iPtBin]);
335 residualList->Add(fHistEPDcaZRes[iPart][iPtBin]);
337 } // loop over pt bins
338 } // loop over particles (pos, neg)
345 //________________________________________________________________________
346 void AliHFEdca::CreateHistogramsPull(TList *pullList){
350 const Int_t nBins = 1000;
351 const Float_t maxXYBin = 20.;
352 const Float_t maxZBin = 20.;
355 // for pull -----------------------------------------------------------------------
356 // fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
357 // fHistDcaZPull[kNParticles][kNPtBins]=0x0;
358 // fHistEPDcaXYPull[kNParticles-2][kNPtBins]=0x0;
359 // fHistEPDcaZPull[kNParticles-2][kNPtBins]=0x0;
362 for(Int_t k=0; k<kNDcaVar; k++){
363 TString histTitle((const char*)fgkPullDcaVarTitle[k]);
365 for(Int_t j=0; j<kNParticles; j++){
366 for(Int_t i=0; i<kNPtBins; i++){
368 TString histName((const char*)fgkParticles[j]);
369 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
370 TString histEPName((const char*)fgkParticles[j]);
371 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
374 fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
375 fHistDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
376 if(j<(kNParticles-2)) {
377 fHistEPDcaXYPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
378 fHistEPDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
381 fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
382 fHistDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
383 if(j<(kNParticles-2)) {
384 fHistEPDcaZPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
385 fHistEPDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
391 // TList *fPullList = 0;
392 pullList->SetOwner();
393 pullList->SetName("pull");
394 for(Int_t iPart=0; iPart<kNParticles; iPart++){
395 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
396 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
397 pullList->Add(fHistDcaZPull[iPart][iPtBin]);
398 if(iPart<(kNParticles-2)){
399 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
400 pullList->Add(fHistDcaZPull[iPart][iPtBin]); }
401 } // loop over pt bins
402 } // loop over particles (pos, neg)
407 //________________________________________________________________________
408 void AliHFEdca::CreateHistogramsDca(TList *dcaList){
410 // define histograms: MC dca
415 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
416 fStat->SetMarkerStyle(20);
417 fStat->SetMarkerColor(3);
418 fStat->SetMarkerSize(1);
421 // fHistDcaXY[kNParticles][kNPtBins]=0x0;
422 // fHistDcaZ[kNParticles][kNPtBins]=0x0;
423 // fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
424 // fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
426 const Int_t nBins = 1000;
427 const Float_t maxXYBin = 1000.;
428 const Float_t maxZBin = 1000.;
431 for(Int_t k=0; k<kNDcaVar; k++){
432 TString histTitle((const char*)fgkDcaVarTitle[k]);
434 for(Int_t j=0; j<kNParticles; j++){
435 for(Int_t i=0; i<kNPtBins; i++){
437 TString histName((const char*)fgkParticles[j]);
438 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
440 TString histNameEP((const char*)fgkParticles[j]);
441 histNameEP += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
444 fHistDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
445 fHistDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
447 if(j<(kNParticles-2)){
448 fHistEPDcaXY[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
449 fHistEPDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);}
452 fHistDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
453 fHistDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
454 if(j<(kNParticles-2)){
455 fHistEPDcaZ[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxZBin, maxZBin);
456 fHistEPDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);}
462 // TList *fDcaList = 0;
464 dcaList->SetName("mcDcaDistr");
465 for(Int_t iPart=0; iPart<kNParticles; iPart++){
466 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
467 dcaList->Add(fHistDcaXY[iPart][iPtBin]);
468 dcaList->Add(fHistDcaZ[iPart][iPtBin]);
469 if(iPart<(kNParticles-2)) {
470 dcaList->Add(fHistEPDcaXY[iPart][iPtBin]);
471 dcaList->Add(fHistEPDcaZ[iPart][iPtBin]); }
472 } // loop over pt bins
473 } // loop over particles (pos, neg)
479 //________________________________________________________________________
480 void AliHFEdca::CreateHistogramsKfDca(TList *kfDcaList){
482 // define histograms: MC dca
487 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
488 fStat->SetMarkerStyle(20);
489 fStat->SetMarkerColor(3);
490 fStat->SetMarkerSize(1);
493 // fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
494 // fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
496 const Int_t nBins = 1000;
497 const Float_t maxXYBin = 1000.;
498 const Float_t maxZBin = 1000.;
501 for(Int_t k=0; k<kNDcaVar; k++){
502 TString histTitle((const char*)fgkDcaVarTitle[k]);
504 for(Int_t j=0; j<kNParticles; j++){
505 for(Int_t i=0; i<kNPtBins; i++){
506 TString histNameKF((const char*)fgkParticles[j]);
507 histNameKF += Form("_MCpid_KF%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
510 fHistKFDcaXY[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
511 fHistKFDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
514 fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
515 fHistKFDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
521 kfDcaList->SetOwner();
522 kfDcaList->SetName("mcKfDcaDistr");
523 for(Int_t iPart=0; iPart<kNParticles; iPart++){
524 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
525 kfDcaList->Add(fHistKFDcaXY[iPart][iPtBin]);
526 kfDcaList->Add(fHistKFDcaZ[iPart][iPtBin]);
527 } // loop over pt bins
528 } // loop over particles (pos, neg)
530 kfDcaList->Add(fStat);
535 //________________________________________________________________________
536 void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
538 // define histograms: real Data
542 // fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
543 // fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
544 // fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
545 // fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
547 const Int_t nBins = 1000;
548 const Float_t maxXYBin = 1000.;
549 const Float_t maxZBin = 1000.;
551 for(Int_t k=0; k<kNDcaVar; k++){
552 TString histTitle((const char*)fgkDcaVarTitle[k]);
553 for(Int_t j=0; j<kNParticles; j++){
554 for(Int_t i=0; i<kNPtBins; i++){
556 TString histName((const char*)fgkParticles[j]);
557 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
559 TString histNameWo((const char*)fgkParticles[j]);
560 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
563 fHistDataDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
564 fHistDataDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
566 fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
567 fHistDataWoDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
570 fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
571 fHistDataDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
573 fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
574 fHistDataWoDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
580 // TList *fDcaList = 0;
581 dataDcaList->SetOwner();
582 dataDcaList->SetName("dataDcaDistr");
583 for(Int_t iPart=0; iPart<kNParticles; iPart++){
584 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
585 dataDcaList->Add(fHistDataDcaXY[iPart][iPtBin]);
586 dataDcaList->Add(fHistDataDcaZ[iPart][iPtBin]);
588 dataDcaList->Add(fHistDataWoDcaXY[iPart][iPtBin]);
589 dataDcaList->Add(fHistDataWoDcaZ[iPart][iPtBin]);
590 } // loop over pt bins
591 } // loop over particles (pos, neg)
597 //________________________________________________________________________
598 void AliHFEdca::CreateHistogramsDataPull(TList *dataPullList){
602 const Int_t nBins = 1000;
603 const Float_t maxXYBin = 20.;
604 const Float_t maxZBin = 20.;
606 // for pull -----------------------------------------------------------------------
607 // fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
608 // fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
610 // fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
611 // fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
614 for(Int_t k=0; k<kNDcaVar; k++){
615 TString histTitle((const char*)fgkPullDataDcaVarTitle[k]);
617 for(Int_t j=0; j<kNParticles; j++){
618 for(Int_t i=0; i<kNPtBins; i++){
620 TString histName((const char*)fgkParticles[j]);
621 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
623 TString histNameWo((const char*)fgkParticles[j]);
624 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
627 fHistDataDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
628 fHistDataDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
630 fHistDataWoDcaXYPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
631 fHistDataWoDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
634 fHistDataDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
635 fHistDataDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
637 fHistDataWoDcaZPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
638 fHistDataWoDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
644 // TList *fDataPullList = 0;
645 dataPullList->SetOwner();
646 dataPullList->SetName("dataPull");
647 for(Int_t iPart=0; iPart<kNParticles; iPart++){
648 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
649 dataPullList->Add(fHistDataDcaXYPull[iPart][iPtBin]);
650 dataPullList->Add(fHistDataDcaZPull[iPart][iPtBin]);
652 dataPullList->Add(fHistDataWoDcaXYPull[iPart][iPtBin]);
653 dataPullList->Add(fHistDataWoDcaZPull[iPart][iPtBin]);
654 } // loop over pt bins
655 } // loop over particles (pos, neg)
659 //________________________________________________________________________
660 void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
662 // define histograms: vertex
666 // fHistMCvertex[kNVertexVar]=0x0;
667 // fHistESDvertex[kNVertexVar]=0x0;
669 const Int_t nBins = 1000;
670 const Float_t minXBin = -0.2e4;
671 const Float_t maxXBin = 0.2e4;
672 const Float_t minYBin = -0.5e4;
673 const Float_t maxYBin = 0.5e4;
674 const Float_t minZBin = -1.5e5;
675 const Float_t maxZBin = 1.5e5;
677 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
678 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
680 for(Int_t k=0; k<kNVertexVar; k++){
681 TString histTitle((const char*)fgkVertexVarTitle[k]);
682 TString histNameMC((const char*)fgkVertexVar[k]);
683 histNameMC += Form("_MC");
684 TString histNameESD((const char*)fgkVertexVar[k]);
685 histNameESD += Form("_ESD");
687 fHistMCvertex[k] = new TH1F((const char*)histNameMC, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
688 fHistMCvertex[k]->SetLineColor(k+2);
690 fHistESDvertex[k] = new TH1F((const char*)histNameESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
691 fHistESDvertex[k]->SetLineColor(k+2);
694 vertexList->SetOwner();
695 vertexList->SetName("vertexDistr");
697 for(Int_t k=0; k<kNVertexVar; k++){
698 vertexList->Add(fHistMCvertex[k]);
699 vertexList->Add(fHistESDvertex[k]);
706 //________________________________________________________________________
707 void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
709 // define histograms: vertex
713 // fHistDatavertex[kNVertexVar]=0x0;
715 const Int_t nBins = 1000;
716 const Float_t minXBin = -0.2e4;
717 const Float_t maxXBin = 0.2e4;
718 const Float_t minYBin = -0.5e4;
719 const Float_t maxYBin = 0.5e4;
720 const Float_t minZBin = -1.5e5;
721 const Float_t maxZBin = 1.5e5;
723 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
724 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
726 for(Int_t k=0; k<kNVertexVar; k++){
727 TString histTitle((const char*)fgkVertexVarTitle[k]);
728 TString histNameDataESD((const char*)fgkVertexVar[k]);
729 histNameDataESD += Form("_data");
731 fHistDatavertex[k] = new TH1F((const char*)histNameDataESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
732 fHistDatavertex[k]->SetLineColor(k+2);
735 // TList *fVDaraVertexList = 0;
736 dataVertexList->SetOwner();
737 dataVertexList->SetName("dataVertexDistr");
739 for(Int_t k=0; k<kNVertexVar; k++){
740 dataVertexList->Add(fHistDatavertex[k]);
745 //_______________________________________________________________________________________________
746 void AliHFEdca::CreateHistogramsPid(TList *mcPidList){
748 // define histograms which fills combined PID
751 const Char_t *mcOResd[2]={"mcPt", "esdPt"};
753 for(Int_t iPart=0; iPart<kNParticles; iPart++){
754 TString histTitleMc((const char*)fgkParticles[iPart]);
755 TString histTitleEsd((const char*)fgkParticles[iPart]);
756 histTitleMc += Form("_McPid_%s;p_{T} [GeV/c];counts", mcOResd[0]);
757 histTitleEsd += Form("_EsdPid_%s;p_{T} [GeV/c];counts", mcOResd[1]);
759 TString histNameMc((const char*)fgkParticles[iPart]);
760 TString histNameEsd((const char*)fgkParticles[iPart]);
761 histNameMc+=Form("_McPid_%s", mcOResd[0]);
762 histNameEsd+=Form("_EsdPid_%s", mcOResd[1]);
764 fHistMcPid[iPart] = new TH1F(histNameMc, histTitleMc, kNPtBins, fgkPtIntv);
765 fHistEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
769 mcPidList->SetOwner();
770 mcPidList->SetName("combinedPid");
772 for(Int_t iPart=0; iPart<kNParticles; iPart++){
773 mcPidList->Add(fHistMcPid[iPart]);
774 mcPidList->Add(fHistEsdPid[iPart]);
779 //_______________________________________________________________________________________________
780 void AliHFEdca::CreateHistogramsDataPid(TList *pidList){
782 // define histograms which fills combined PID: data
787 for(Int_t iPart=0; iPart<kNParticles; iPart++){
788 TString histTitleEsd((const char*)fgkParticles[iPart]);
789 histTitleEsd+=Form("_DataEsdPid_esdPt;p_{T} [GeV/c];counts");
790 TString histNameEsd((const char*)fgkParticles[iPart]);
791 histNameEsd+=Form("_DataEsdPid");
793 fHistDataEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
798 pidList->SetName("dataCombinedPid");
800 for(Int_t iPart=0; iPart<kNParticles; iPart++)
801 pidList->Add(fHistDataEsdPid[iPart]);
808 //_______________________________________________________________________________________________
809 void AliHFEdca::FillHistogramsDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
813 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
815 mcPrimV[0] = mcPrimVtx->GetX();
816 mcPrimV[1] = mcPrimVtx->GetY();
817 mcPrimV[2] = mcPrimVtx->GetZ();
819 Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
821 // filling historgams track by track
822 // obtaining reconstructed dca ------------------------------------------------------------------
824 Float_t esdpx = track->Px();
825 Float_t esdpy = track->Py();
826 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
828 // obtaining errors of dca ------------------------------------------------------------------
829 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
831 primV[0] = primVtx->GetXv();
832 primV[1] = primVtx->GetYv();
833 primV[2] = primVtx->GetZv();
835 Float_t magneticField = 0; // initialized as 5kG
836 magneticField = esdEvent->GetMagneticField(); // in kG
838 Double_t beampiperadius=3.;
839 Double_t dz[2]; // error of dca in cm
842 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
844 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
846 TParticle *part = mctrack->Particle();
848 Float_t vx = part->Vx(); // in cm
849 Float_t vy = part->Vy(); // in cm
850 Float_t vz = part->Vz(); // in cm
852 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
854 Float_t mcpx = part->Px();
855 Float_t mcpy = part->Py();
856 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
858 Int_t pdg = part->GetPdgCode();
859 Int_t esdPid = GetCombinedPid(track);
862 if(pdg==kPDGelectron || pdg==kPDGmuon
863 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
865 // calculate mcDca ------------------------------------------------------------------
866 const Float_t conv[2] = {1.783/1.6, 2.99792458};
867 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
869 Float_t nx = esdpx/mcpt;
870 Float_t ny = esdpy/mcpt;
873 radius = TMath::Abs(radiusMc);
875 Double_t dxy = vxy - mcVtxXY; // in cm
876 Double_t dvx = vx - mcPrimV[0]; // in cm
877 Double_t dvy = vy - mcPrimV[1]; // in cm
879 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
881 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
882 Double_t residual[2] = {0, 0};
883 Double_t pull[2] = {0, 0};
884 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
885 for(Int_t i=0; i<2; i++){
886 residual[i] = dz[i] - mcDca[i]; // in centimeters
887 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
891 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
893 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
894 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
895 if(pdg==fgkPdgParticle[iPart]) {
896 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
897 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
898 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
899 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
900 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
901 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
904 if(esdPid==fgkPdgParticle[iPart]) {
905 fHistEPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
906 fHistEPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
907 fHistEPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
908 fHistEPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
909 fHistEPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
910 fHistEPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
918 } // particle id loop
920 // for charged particles: no pid
921 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
922 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
924 if(charge>0) iPart = 11;
925 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
926 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
927 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
928 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
929 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
930 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
938 //_______________________________________________________________________________________________
939 void AliHFEdca::FillHistogramsKfDca(AliESDEvent * const esdEvent, AliESDtrack * const track, const AliMCEvent * const mcEvent)
943 // filling historgams track by track
945 // obtaining reconstructed dca ------------------------------------------------------------------
946 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
947 Float_t magneticField = 0; // initialized as 5kG
948 magneticField = esdEvent->GetMagneticField(); // in kG
950 Float_t esdpx = track->Px();
951 Float_t esdpy = track->Py();
952 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
954 Int_t charge = track->Charge();
956 Double_t beampiperadius=3.;
957 Double_t dz[2]; // error of dca in cm
959 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
961 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
963 TParticle *part = mctrack->Particle();
964 Int_t pdg = part->GetPdgCode();
966 // calculate dca using AliKFParticle class------------------------------------------------------------------
967 Double_t kfdz[3] = {0, 0, 0};
968 Double_t kfdzwith[3] = {0, 0, 0};
970 Int_t trkID = track->GetID();
972 AliKFParticle::SetField(magneticField);
973 AliKFParticle kfParticle(*track, pdg);
975 // prepare kfprimary vertex
976 AliKFVertex kfESDprimary;
977 // Reconstruct Primary Vertex (with ESD tracks)
978 Int_t n=primVtx->GetNIndices();
979 if (n>0 && primVtx->GetStatus()){
980 kfESDprimary = AliKFVertex(*primVtx);
982 Double_t dcaXYWithTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
983 Double_t dcaWithTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
984 Double_t dcaZWithTrk = 0;
985 if(TMath::Abs(dcaWithTrk)>=TMath::Abs(dcaXYWithTrk))
986 dcaZWithTrk =TMath::Sqrt(dcaWithTrk*dcaWithTrk-dcaXYWithTrk*dcaXYWithTrk)*((dz[1]*-1<=0)?1:-1);
987 kfdzwith[0] = dcaXYWithTrk;
988 kfdzwith[1] = dcaZWithTrk;
989 kfdzwith[2] = dcaWithTrk; // with current track
991 Double_t dcaXYWoTrk = 0;
992 Double_t dcaZWoTrk = 0;
993 Double_t dcaWoTrk = 0;
995 UShort_t *priIndex = primVtx->GetIndices();
997 for (Int_t i=0;i<n;i++){
999 Int_t idx = Int_t(priIndex[i]);
1001 kfESDprimary -= kfParticle;
1002 dcaXYWoTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
1003 dcaWoTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
1004 if((dcaWoTrk-dcaXYWoTrk)>=0)
1005 dcaZWoTrk = TMath::Abs(dcaWoTrk*dcaWoTrk - dcaXYWoTrk*dcaXYWoTrk)*((dz[1]*-1<=0)?1:-1);
1006 } // remove current track from this calculation
1007 } // loop over all primary vertex contributors
1010 kfdz[0] = dcaXYWoTrk;
1011 kfdz[1] = dcaZWoTrk;
1014 } // only if n contributor > 0 and primVtx constructed
1018 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
1019 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
1020 if(kfdzwith[0]==0 && dz[0]!=0) fStat->Fill(3); // 0 from KF particle (with current track)
1022 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
1023 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
1024 if(kfdz[0]==0 && dz[0]!=0) fStat->Fill(-3); // 0 from KF particle (without current track)
1026 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1028 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1029 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1030 if(pdg==fgkPdgParticle[iPart]) {
1031 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1.0e4);
1032 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1.0e4);
1039 } // particle id loop
1041 // for charged particles: no pid
1042 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1043 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1045 if(charge>0) iPart = 11;
1046 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1047 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1057 //_______________________________________________________________________________________________
1058 void AliHFEdca::FillHistogramsVtx(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
1062 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1063 Double_t mcPrimV[3];
1064 mcPrimV[0] = mcPrimVtx->GetX();
1065 mcPrimV[1] = mcPrimVtx->GetY();
1066 mcPrimV[2] = mcPrimVtx->GetZ();
1068 // obtaining errors of dca ------------------------------------------------------------------
1069 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1071 primV[0] = primVtx->GetXv();
1072 primV[1] = primVtx->GetYv();
1073 primV[2] = primVtx->GetZv();
1075 for(Int_t i=0; i<kNVertexVar; i++){
1076 fHistMCvertex[i]->Fill(mcPrimV[i]*1.0e4);
1077 fHistESDvertex[i]->Fill(primV[i]*1.0e4);
1082 //_______________________________________________________________________________________________
1083 void AliHFEdca::FillHistogramsPid(AliESDtrack * const track, const AliMCEvent * const mcEvent)
1087 // filling historgams track by track
1088 Float_t esdpx = track->Px();
1089 Float_t esdpy = track->Py();
1090 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1092 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1093 if(!mctrack) return;
1094 TParticle *part = mctrack->Particle();
1096 Float_t mcpx = part->Px();
1097 Float_t mcpy = part->Py();
1098 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1100 Int_t pdg = part->GetPdgCode();
1101 Int_t esdPid = GetCombinedPid(track);
1104 Double_t ptMom[2] = {mcpt, esdpt};
1106 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1107 if(pdg==fgkPdgParticle[iPart]) // pid all by MC
1108 fHistMcPid[iPart]->Fill(ptMom[0]);
1110 if(esdPid==fgkPdgParticle[iPart]) // pid all by combined pid
1111 fHistEsdPid[iPart]->Fill(ptMom[1]);
1112 } // loop over particles
1115 if(pdg==kPDGelectron || pdg==kPDGmuon || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton)
1116 fHistMcPid[10]->Fill(ptMom[0]);
1117 if(pdg==-kPDGelectron || pdg==-kPDGmuon || pdg==kPDGpion || pdg==kPDGkaon || pdg==kPDGproton)
1118 fHistMcPid[11]->Fill(ptMom[0]);
1119 if(esdPid==kPDGelectron || esdPid==kPDGmuon || esdPid==-kPDGpion || esdPid==-kPDGkaon || esdPid==-kPDGproton)
1120 fHistEsdPid[10]->Fill(ptMom[1]);
1121 if(esdPid==-kPDGelectron || esdPid==-kPDGmuon || esdPid==kPDGpion || esdPid==kPDGkaon || esdPid==kPDGproton)
1122 fHistEsdPid[11]->Fill(ptMom[1]);
1126 ////_______________________________________________________________________________________________
1127 void AliHFEdca::FillHistogramsDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1129 // filling historgams track by track
1130 // obtaining reconstructed dca --------------------------------------------------------------
1132 Float_t esdpx = track->Px();
1133 Float_t esdpy = track->Py();
1134 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1135 Int_t charge = track->Charge();
1137 // obtaining errors of dca ------------------------------------------------------------------
1138 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1140 primV[0] = primVtx->GetXv();
1141 primV[1] = primVtx->GetYv();
1142 primV[2] = primVtx->GetZv();
1145 Float_t magneticField = 0; // initialized as 5kG
1146 magneticField = esdEvent->GetMagneticField(); // in kG
1148 Double_t beampiperadius=3.;
1149 Double_t dz[2]; // error of dca in cm
1150 Double_t covardz[3];
1152 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1155 Double_t pull[2] = {0, 0};
1156 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1157 for(Int_t i=0; i<2; i++){
1158 if(error[i]!=0)pull[i] = dz[i]/error[i]; // unitless
1161 // get dca when current track is not included
1163 Double_t dzwo[2], covardzwo[3];
1164 Double_t pullwo[2] = {0, 0};
1165 if(!track->PropagateToDCA(vtxESDSkip, magneticField, beampiperadius, dzwo, covardzwo)) return; // protection
1167 Double_t errorwo[2] ={TMath::Sqrt(TMath::Abs(covardzwo[0])), TMath::Sqrt(TMath::Abs(covardzwo[2]))};
1168 for(Int_t i=0; i<2; i++){
1169 if(errorwo[i]!=0) pullwo[i] = dzwo[i]/errorwo[i]; // unitless
1173 Int_t esdPid = GetCombinedPid(track);
1175 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1177 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1178 if(esdPid==fgkPdgParticle[iPart] && (esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1])) {
1179 fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1180 fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1181 fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1182 fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1183 // w/o current track
1184 fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1185 fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1186 fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1187 fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1194 // for charged particles
1195 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1196 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1198 if(charge>0) iPart = 11;
1199 fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1200 fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1201 fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1202 fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1203 // without current track
1204 fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1205 fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1206 fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1207 fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1216 //_______________________________________________________________________________________________
1217 void AliHFEdca::FillHistogramsDataVtx(AliESDEvent * const esdEvent)
1221 // obtaining errors of dca ------------------------------------------------------------------
1222 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1224 primV[0] = primVtx->GetXv();
1225 primV[1] = primVtx->GetYv();
1226 primV[2] = primVtx->GetZv();
1228 // require events with at least 3 contributors for primary vertex construction
1229 Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1230 if(nEsdPrimVtxCtrb<1) return; // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1231 for(Int_t i=0; i<kNVertexVar; i++)
1232 fHistDatavertex[i]->Fill(primV[i]*1e4);
1237 ////_______________________________________________________________________________________________
1238 void AliHFEdca::FillHistogramsDataPid(AliESDtrack * const track)
1240 // filling historgams track by track
1241 // obtaining reconstructed dca --------------------------------------------------------------
1243 Float_t esdpx = track->Px();
1244 Float_t esdpy = track->Py();
1245 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1246 Int_t charge = track->Charge();
1248 Int_t esdPid = GetCombinedPid(track);
1250 for(Int_t iPart=0; iPart<kNParticles; iPart++){
1251 if(iPart<(kNParticles-2)){
1252 if(esdPid==fgkPdgParticle[iPart]) fHistDataEsdPid[iPart]->Fill(esdpt);
1256 if(charge<0) fHistDataEsdPid[10]->Fill(esdpt);
1257 if(charge>0) fHistDataEsdPid[11]->Fill(esdpt);
1262 //_________________________________________________________________________________________________
1263 void AliHFEdca::ApplyExtraCuts(AliESDEvent * const esdEvent, Int_t nMinPrimVtxContributor)
1267 // only one extra cut, number of contributors to each primary vertex
1270 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1271 Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1272 if(nEsdPrimVtxCtrb<nMinPrimVtxContributor) return;
1273 // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1277 //_____________________________________________________
1278 Int_t AliHFEdca::GetCombinedPid(AliESDtrack *const track)
1281 // combined detector pid
1282 Double_t prob[AliPID::kSPECIES];
1283 track->GetESDpid(prob);
1286 Double_t priors[AliPID::kSPECIESN];
1293 Int_t charge = track->Charge();
1297 pid.SetPriors(priors);
1298 pid.SetProbabilities(prob);
1300 // identify particle as the most probable
1302 Double_t pelectron = pid.GetProbability(AliPID::kElectron);
1303 if(pelectron > pid.GetProbability(AliPID::kMuon) &&
1304 pelectron > pid.GetProbability(AliPID::kPion) &&
1305 pelectron > pid.GetProbability(AliPID::kKaon) &&
1306 pelectron > pid.GetProbability(AliPID::kProton) ) esdPid = -kPDGelectron;
1308 Double_t pmuon = pid.GetProbability(AliPID::kMuon);
1309 if(pmuon > pid.GetProbability(AliPID::kElectron) &&
1310 pmuon > pid.GetProbability(AliPID::kPion) &&
1311 pmuon > pid.GetProbability(AliPID::kKaon) &&
1312 pmuon > pid.GetProbability(AliPID::kProton) ) esdPid = -kPDGmuon;
1314 Double_t ppion = pid.GetProbability(AliPID::kPion);
1315 if(ppion > pid.GetProbability(AliPID::kElectron) &&
1316 ppion > pid.GetProbability(AliPID::kMuon) &&
1317 ppion > pid.GetProbability(AliPID::kKaon) &&
1318 ppion > pid.GetProbability(AliPID::kProton) ) esdPid = kPDGpion;
1320 Double_t pkaon = pid.GetProbability(AliPID::kKaon);
1321 if(pkaon > pid.GetProbability(AliPID::kElectron) &&
1322 pkaon > pid.GetProbability(AliPID::kMuon) &&
1323 pkaon > pid.GetProbability(AliPID::kPion) &&
1324 pkaon > pid.GetProbability(AliPID::kProton) ) esdPid = kPDGkaon;
1326 Double_t pproton = pid.GetProbability(AliPID::kProton);
1327 if(pproton > pid.GetProbability(AliPID::kElectron) &&
1328 pproton > pid.GetProbability(AliPID::kMuon) &&
1329 pproton > pid.GetProbability(AliPID::kPion) &&
1330 pproton > pid.GetProbability(AliPID::kKaon) ) esdPid = kPDGproton;
1333 return charge*esdPid;
1340 //________________________________________________________________________
1341 void AliHFEdca::CreateHistogramsHfeDca(TList *hfeDcaList){
1343 // define histograms: hfe pid electrons in MC
1346 const Int_t nBinsDca = 1000;
1347 const Float_t maxXYBinDca = 1000.;
1348 const Float_t maxZBinDca = 1000.;
1350 const Int_t nBinsPull = 1000;
1351 const Float_t maxXYBinPull = 20.;
1352 const Float_t maxZBinPull = 20.;
1354 const Char_t *mcOResd[2]={"mcPt", "esdPt"};
1356 // fHistHPDcaXY[2][kNPtBins]=0x0;
1357 // fHistHPDcaZ[2][kNPtBins]=0x0;
1358 // fHistHPDcaXYRes[2][kNPtBins]=0x0;
1359 // fHistHPDcaZRes[2][kNPtBins]=0x0;
1360 // fHistHPDcaXYPull[2][kNPtBins]=0x0;
1361 // fHistHPDcaZPull[2][kNPtBins]=0x0;
1364 for(Int_t k=0; k<kNDcaVar; k++){
1366 TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1367 TString histTitleRes((const char*)fgkPullDcaVarTitle[k]);
1368 TString histTitlePull((const char*)fgkResDcaVarTitle[k]);
1370 for(Int_t iPart=0; iPart<2; iPart++){
1371 for(Int_t i=0; i<kNPtBins; i++){
1372 TString histHPName((const char*)fgkParticles[iPart*5]);
1373 histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1374 TString histHPNameRes((const char*)fgkParticles[iPart*5]);
1375 histHPNameRes += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1376 TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1377 histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1380 fHistHPDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1381 fHistHPDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1382 fHistHPDcaXYRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1383 fHistHPDcaXYRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1384 fHistHPDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1385 fHistHPDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1389 fHistHPDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1390 fHistHPDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1391 fHistHPDcaZRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1392 fHistHPDcaZRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1393 fHistHPDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxZBinPull, 1+maxZBinPull);
1394 fHistHPDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1400 // fHistHfePid[2][2] = 0x0; //! HFE pid
1401 for(Int_t id=0; id<2; id++){
1402 for(Int_t iPart=0; iPart<2; iPart++){
1403 TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1404 histTitleHfe+=Form("_MC_HfePid_esdPt;p_{T} [GeV/c];counts");
1405 TString histNameHfe((const char*)fgkParticles[iPart*5]);
1406 histNameHfe+=Form("_MC_HfePid_%s", mcOResd[id]);
1407 fHistHfePid[id][iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1411 hfeDcaList->SetOwner();
1412 hfeDcaList->SetName("hfeDca");
1413 for(Int_t iPart=0; iPart<2; iPart++){
1414 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1415 hfeDcaList->Add(fHistHPDcaXY[iPart][iPtBin]);
1416 hfeDcaList->Add(fHistHPDcaZ[iPart][iPtBin]);
1417 hfeDcaList->Add(fHistHPDcaXYRes[iPart][iPtBin]);
1418 hfeDcaList->Add(fHistHPDcaZRes[iPart][iPtBin]);
1419 hfeDcaList->Add(fHistHPDcaXYPull[iPart][iPtBin]);
1420 hfeDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1422 for(Int_t id=0; id<2; id++)
1423 hfeDcaList->Add(fHistHfePid[id][iPart]);
1429 //________________________________________________________________________
1430 void AliHFEdca::CreateHistogramsHfeDataDca(TList *hfeDataDcaList){
1432 // define histograms: hfe pid electrons in data
1435 const Int_t nBinsDca = 1000;
1436 const Float_t maxXYBinDca = 1000.;
1437 const Float_t maxZBinDca = 1000.;
1439 const Int_t nBinsPull = 1000;
1440 const Float_t maxXYBinPull = 20.;
1441 const Float_t maxZBinPull = 20.;
1444 // fHistHPDataDcaXY[2][kNPtBins]=0x0;
1445 // fHistHPDataDcaZ[2][kNPtBins]=0x0;
1446 // fHistHPDataDcaXYPull[2][kNPtBins]=0x0;
1447 // fHistHPDataDcaZPull[2][kNPtBins]=0x0;
1449 for(Int_t k=0; k<kNDcaVar; k++){
1450 TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1451 TString histTitlePull((const char*)fgkPullDcaVarTitle[k]);
1452 for(Int_t iPart=0; iPart<2; iPart++){
1453 for(Int_t i=0; i<kNPtBins; i++){
1454 TString histHPName((const char*)fgkParticles[iPart*5]);
1455 histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1456 TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1457 histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1460 fHistHPDataDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1461 fHistHPDataDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1462 fHistHPDataDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1463 fHistHPDataDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1467 fHistHPDataDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1468 fHistHPDataDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1469 fHistHPDataDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsDca, 1-maxZBinPull, 1+maxZBinPull);
1470 fHistHPDataDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
1474 } // 2 particle type
1477 //fHistDataHfePid[2] = 0x0; //! HFE pid
1478 for(Int_t iPart=0; iPart<2; iPart++){
1479 TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1480 histTitleHfe+=Form("_Data_HfePid_esdPt;p_{T} [GeV/c];counts");
1481 TString histNameHfe((const char*)fgkParticles[iPart*5]);
1482 histNameHfe+=Form("_Data_HfePid");
1483 fHistDataHfePid[iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1487 hfeDataDcaList->SetOwner();
1488 hfeDataDcaList->SetName("hfeDataDca");
1489 for(Int_t iPart=0; iPart<2; iPart++){
1490 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1491 hfeDataDcaList->Add(fHistHPDataDcaXY[iPart][iPtBin]);
1492 hfeDataDcaList->Add(fHistHPDataDcaZ[iPart][iPtBin]);
1493 hfeDataDcaList->Add(fHistHPDataDcaXYPull[iPart][iPtBin]);
1494 hfeDataDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1496 hfeDataDcaList->Add(fHistDataHfePid[iPart]);
1504 //_______________________________________________________________________________________________
1505 void AliHFEdca::FillHistogramsHfeDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
1507 // the kHFEpid plugin
1509 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1510 Double_t mcPrimV[3];
1511 mcPrimV[0] = mcPrimVtx->GetX();
1512 mcPrimV[1] = mcPrimVtx->GetY();
1513 mcPrimV[2] = mcPrimVtx->GetZ();
1515 Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
1517 // filling historgams track by track
1518 // obtaining reconstructed dca ------------------------------------------------------------------
1519 Float_t esdpx = track->Px();
1520 Float_t esdpy = track->Py();
1521 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1523 // obtaining errors of dca ------------------------------------------------------------------
1524 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1526 primV[0] = primVtx->GetXv();
1527 primV[1] = primVtx->GetYv();
1528 primV[2] = primVtx->GetZv();
1530 Float_t magneticField = 0; // initialized as 5kG
1531 magneticField = esdEvent->GetMagneticField(); // in kG
1533 Double_t beampiperadius=3.;
1534 Double_t dz[2]; // error of dca in cm
1535 Double_t covardz[3];
1536 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1538 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1539 if(!mctrack) return;
1540 TParticle *part = mctrack->Particle();
1542 Float_t vx = part->Vx(); // in cm
1543 Float_t vy = part->Vy(); // in cm
1544 Float_t vz = part->Vz(); // in cm
1546 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1548 Float_t mcpx = part->Px();
1549 Float_t mcpy = part->Py();
1550 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1552 Int_t pdg = part->GetPdgCode();
1555 if(pdg==kPDGelectron || pdg==kPDGmuon
1556 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1558 // calculate mcDca ------------------------------------------------------------------
1559 const Float_t conv[2] = {1.783/1.6, 2.99792458};
1560 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
1562 Float_t nx = esdpx/mcpt;
1563 Float_t ny = esdpy/mcpt;
1566 radius = TMath::Abs(radiusMc);
1568 Double_t dxy = vxy - mcVtxXY; // in cm
1569 Double_t dvx = vx - mcPrimV[0]; // in cm
1570 Double_t dvy = vy - mcPrimV[1]; // in cm
1572 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
1574 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
1575 Double_t residual[2] = {0, 0};
1576 Double_t pull[2] = {0, 0};
1577 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1578 for(Int_t i=0; i<2; i++){
1579 residual[i] = dz[i] - mcDca[i]; // in centimeters
1580 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
1584 if(track->Charge()<0) iPart = 0; // electron
1585 if(track->Charge()>0) iPart = 1; // positron
1586 if(track->Charge()==0) {
1587 printf("this is not an electron! Check HFEpid method");
1590 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1591 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1592 fHistHPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
1593 fHistHPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
1594 fHistHPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1595 fHistHPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1596 fHistHPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
1597 fHistHPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
1605 fHistHfePid[iPart][0]->Fill(esdpt);
1606 fHistHfePid[iPart][1]->Fill(mcpt);
1611 //_______________________________________________________________________________________________
1612 void AliHFEdca::FillHistogramsHfeDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1614 // filling historgams track by track
1615 // obtaining reconstructed dca --------------------------------------------------------------
1617 Float_t esdpx = track->Px();
1618 Float_t esdpy = track->Py();
1619 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1620 Int_t charge = track->Charge();
1622 // obtaining errors of dca ------------------------------------------------------------------
1623 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!
1625 primV[0] = primVtx->GetXv();
1626 primV[1] = primVtx->GetYv();
1627 primV[2] = primVtx->GetZv();
1629 Float_t magneticField = 0; // initialized as 5kG
1630 magneticField = esdEvent->GetMagneticField(); // in kG
1631 Double_t beampiperadius=3.;
1633 Double_t dz[2]; // error of dca in cm
1634 Double_t covardz[3];
1636 if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
1638 Double_t pull[2] = {0, 0};
1639 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1640 for(Int_t i=0; i<2; i++){
1641 if(error[i]!=0) pull[i] = dz[i]/error[i]; // unitless
1645 if(charge<0) iPart = 0; // electron
1646 if(charge>0) iPart = 1; // positron
1648 printf("this is not an electron! Check HFEpid method\n");
1652 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1653 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]) {
1654 fHistHPDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1655 fHistHPDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1656 fHistHPDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1657 fHistHPDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1663 fHistDataHfePid[iPart]->Fill(esdpt);