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()const{
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(const AliESDEvent * const esdEvent, const 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 AliESDtrack ctrack(*track);
843 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
845 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
847 TParticle *part = mctrack->Particle();
849 Float_t vx = part->Vx(); // in cm
850 Float_t vy = part->Vy(); // in cm
851 Float_t vz = part->Vz(); // in cm
853 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
855 Float_t mcpx = part->Px();
856 Float_t mcpy = part->Py();
857 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
859 Int_t pdg = part->GetPdgCode();
860 Int_t esdPid = GetCombinedPid(track);
863 if(pdg==kPDGelectron || pdg==kPDGmuon
864 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
866 // calculate mcDca ------------------------------------------------------------------
867 const Float_t conv[2] = {1.783/1.6, 2.99792458};
868 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
870 Float_t nx = esdpx/mcpt;
871 Float_t ny = esdpy/mcpt;
874 radius = TMath::Abs(radiusMc);
876 Double_t dxy = vxy - mcVtxXY; // in cm
877 Double_t dvx = vx - mcPrimV[0]; // in cm
878 Double_t dvy = vy - mcPrimV[1]; // in cm
880 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
882 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
883 Double_t residual[2] = {0, 0};
884 Double_t pull[2] = {0, 0};
885 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
886 for(Int_t i=0; i<2; i++){
887 residual[i] = dz[i] - mcDca[i]; // in centimeters
888 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
892 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
894 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
895 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
896 if(pdg==fgkPdgParticle[iPart]) {
897 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
898 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
899 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
900 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
901 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
902 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
905 if(esdPid==fgkPdgParticle[iPart]) {
906 fHistEPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
907 fHistEPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
908 fHistEPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
909 fHistEPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
910 fHistEPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
911 fHistEPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
919 } // particle id loop
921 // for charged particles: no pid
922 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
923 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
925 if(charge>0) iPart = 11;
926 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
927 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
928 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
929 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
930 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
931 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
939 //_______________________________________________________________________________________________
940 void AliHFEdca::FillHistogramsKfDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliMCEvent * const mcEvent)
944 // filling historgams track by track
946 // obtaining reconstructed dca ------------------------------------------------------------------
947 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
948 Float_t magneticField = 0; // initialized as 5kG
949 magneticField = esdEvent->GetMagneticField(); // in kG
951 Float_t esdpx = track->Px();
952 Float_t esdpy = track->Py();
953 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
955 Int_t charge = track->Charge();
957 Double_t beampiperadius=3.;
958 Double_t dz[2]; // error of dca in cm
960 AliESDtrack ctrack(*track);
961 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
963 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
965 TParticle *part = mctrack->Particle();
966 Int_t pdg = part->GetPdgCode();
968 // calculate dca using AliKFParticle class------------------------------------------------------------------
969 Double_t kfdz[3] = {0, 0, 0};
970 Double_t kfdzwith[3] = {0, 0, 0};
972 Int_t trkID = track->GetID();
974 AliKFParticle::SetField(magneticField);
975 AliKFParticle kfParticle(*track, pdg);
977 // prepare kfprimary vertex
978 AliKFVertex kfESDprimary;
979 // Reconstruct Primary Vertex (with ESD tracks)
980 Int_t n=primVtx->GetNIndices();
981 if (n>0 && primVtx->GetStatus()){
982 kfESDprimary = AliKFVertex(*primVtx);
984 Double_t dcaXYWithTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
985 Double_t dcaWithTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
986 Double_t dcaZWithTrk = 0;
987 if(TMath::Abs(dcaWithTrk)>=TMath::Abs(dcaXYWithTrk))
988 dcaZWithTrk =TMath::Sqrt(dcaWithTrk*dcaWithTrk-dcaXYWithTrk*dcaXYWithTrk)*((dz[1]*-1<=0)?1:-1);
989 kfdzwith[0] = dcaXYWithTrk;
990 kfdzwith[1] = dcaZWithTrk;
991 kfdzwith[2] = dcaWithTrk; // with current track
993 Double_t dcaXYWoTrk = 0;
994 Double_t dcaZWoTrk = 0;
995 Double_t dcaWoTrk = 0;
997 UShort_t *priIndex = primVtx->GetIndices();
999 for (Int_t i=0;i<n;i++){
1001 Int_t idx = Int_t(priIndex[i]);
1003 kfESDprimary -= kfParticle;
1004 dcaXYWoTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
1005 dcaWoTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
1006 if((dcaWoTrk-dcaXYWoTrk)>=0)
1007 dcaZWoTrk = TMath::Abs(dcaWoTrk*dcaWoTrk - dcaXYWoTrk*dcaXYWoTrk)*((dz[1]*-1<=0)?1:-1);
1008 } // remove current track from this calculation
1009 } // loop over all primary vertex contributors
1012 kfdz[0] = dcaXYWoTrk;
1013 kfdz[1] = dcaZWoTrk;
1016 } // only if n contributor > 0 and primVtx constructed
1020 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
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(2); // swapped sign
1022 if(kfdzwith[0]==0 && dz[0]!=0) fStat->Fill(3); // 0 from KF particle (with current track)
1024 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
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(-2); // swapped sign
1026 if(kfdz[0]==0 && dz[0]!=0) fStat->Fill(-3); // 0 from KF particle (without current track)
1028 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1030 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1031 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1032 if(pdg==fgkPdgParticle[iPart]) {
1033 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1.0e4);
1034 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1.0e4);
1041 } // particle id loop
1043 // for charged particles: no pid
1044 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1045 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1047 if(charge>0) iPart = 11;
1048 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1049 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1059 //_______________________________________________________________________________________________
1060 void AliHFEdca::FillHistogramsVtx(const AliESDEvent *const esdEvent, const AliMCEvent *const mcEvent)
1064 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1065 Double_t mcPrimV[3];
1066 mcPrimV[0] = mcPrimVtx->GetX();
1067 mcPrimV[1] = mcPrimVtx->GetY();
1068 mcPrimV[2] = mcPrimVtx->GetZ();
1070 // obtaining errors of dca ------------------------------------------------------------------
1071 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1073 primV[0] = primVtx->GetXv();
1074 primV[1] = primVtx->GetYv();
1075 primV[2] = primVtx->GetZv();
1077 for(Int_t i=0; i<kNVertexVar; i++){
1078 fHistMCvertex[i]->Fill(mcPrimV[i]*1.0e4);
1079 fHistESDvertex[i]->Fill(primV[i]*1.0e4);
1084 //_______________________________________________________________________________________________
1085 void AliHFEdca::FillHistogramsPid(const AliESDtrack * const track, const AliMCEvent * const mcEvent)
1089 // filling historgams track by track
1090 Float_t esdpx = track->Px();
1091 Float_t esdpy = track->Py();
1092 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1094 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1095 if(!mctrack) return;
1096 TParticle *part = mctrack->Particle();
1098 Float_t mcpx = part->Px();
1099 Float_t mcpy = part->Py();
1100 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1102 Int_t pdg = part->GetPdgCode();
1103 Int_t esdPid = GetCombinedPid(track);
1106 Double_t ptMom[2] = {mcpt, esdpt};
1108 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1109 if(pdg==fgkPdgParticle[iPart]) // pid all by MC
1110 fHistMcPid[iPart]->Fill(ptMom[0]);
1112 if(esdPid==fgkPdgParticle[iPart]) // pid all by combined pid
1113 fHistEsdPid[iPart]->Fill(ptMom[1]);
1114 } // loop over particles
1117 if(pdg==kPDGelectron || pdg==kPDGmuon || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton)
1118 fHistMcPid[10]->Fill(ptMom[0]);
1119 if(pdg==-kPDGelectron || pdg==-kPDGmuon || pdg==kPDGpion || pdg==kPDGkaon || pdg==kPDGproton)
1120 fHistMcPid[11]->Fill(ptMom[0]);
1121 if(esdPid==kPDGelectron || esdPid==kPDGmuon || esdPid==-kPDGpion || esdPid==-kPDGkaon || esdPid==-kPDGproton)
1122 fHistEsdPid[10]->Fill(ptMom[1]);
1123 if(esdPid==-kPDGelectron || esdPid==-kPDGmuon || esdPid==kPDGpion || esdPid==kPDGkaon || esdPid==kPDGproton)
1124 fHistEsdPid[11]->Fill(ptMom[1]);
1128 ////_______________________________________________________________________________________________
1129 void AliHFEdca::FillHistogramsDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
1131 // filling historgams track by track
1132 // obtaining reconstructed dca --------------------------------------------------------------
1134 Float_t esdpx = track->Px();
1135 Float_t esdpy = track->Py();
1136 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1137 Int_t charge = track->Charge();
1139 // obtaining errors of dca ------------------------------------------------------------------
1140 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1142 primV[0] = primVtx->GetXv();
1143 primV[1] = primVtx->GetYv();
1144 primV[2] = primVtx->GetZv();
1147 Float_t magneticField = 0; // initialized as 5kG
1148 magneticField = esdEvent->GetMagneticField(); // in kG
1150 Double_t beampiperadius=3.;
1151 Double_t dz[2]; // error of dca in cm
1152 Double_t covardz[3];
1154 AliESDtrack ctrack(*track);
1155 if(!ctrack.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(!ctrack.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(const 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(const 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(const 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(const 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(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const 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 AliESDtrack ctrack(*track);
1540 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1542 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1543 if(!mctrack) return;
1544 TParticle *part = mctrack->Particle();
1546 Float_t vx = part->Vx(); // in cm
1547 Float_t vy = part->Vy(); // in cm
1548 Float_t vz = part->Vz(); // in cm
1550 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1552 Float_t mcpx = part->Px();
1553 Float_t mcpy = part->Py();
1554 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1556 Int_t pdg = part->GetPdgCode();
1559 if(pdg==kPDGelectron || pdg==kPDGmuon
1560 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1562 // calculate mcDca ------------------------------------------------------------------
1563 const Float_t conv[2] = {1.783/1.6, 2.99792458};
1564 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
1566 Float_t nx = esdpx/mcpt;
1567 Float_t ny = esdpy/mcpt;
1570 radius = TMath::Abs(radiusMc);
1572 Double_t dxy = vxy - mcVtxXY; // in cm
1573 Double_t dvx = vx - mcPrimV[0]; // in cm
1574 Double_t dvy = vy - mcPrimV[1]; // in cm
1576 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
1578 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
1579 Double_t residual[2] = {0, 0};
1580 Double_t pull[2] = {0, 0};
1581 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1582 for(Int_t i=0; i<2; i++){
1583 residual[i] = dz[i] - mcDca[i]; // in centimeters
1584 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
1588 if(track->Charge()<0) iPart = 0; // electron
1589 if(track->Charge()>0) iPart = 1; // positron
1590 if(track->Charge()==0) {
1591 printf("this is not an electron! Check HFEpid method");
1594 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1595 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1596 fHistHPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
1597 fHistHPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
1598 fHistHPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1599 fHistHPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1600 fHistHPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
1601 fHistHPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
1609 fHistHfePid[iPart][0]->Fill(esdpt);
1610 fHistHfePid[iPart][1]->Fill(mcpt);
1615 //_______________________________________________________________________________________________
1616 void AliHFEdca::FillHistogramsHfeDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
1618 // filling historgams track by track
1619 // obtaining reconstructed dca --------------------------------------------------------------
1621 Float_t esdpx = track->Px();
1622 Float_t esdpy = track->Py();
1623 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1624 Int_t charge = track->Charge();
1626 // obtaining errors of dca ------------------------------------------------------------------
1627 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!
1629 primV[0] = primVtx->GetXv();
1630 primV[1] = primVtx->GetYv();
1631 primV[2] = primVtx->GetZv();
1633 Float_t magneticField = 0; // initialized as 5kG
1634 magneticField = esdEvent->GetMagneticField(); // in kG
1635 Double_t beampiperadius=3.;
1637 Double_t dz[2]; // error of dca in cm
1638 Double_t covardz[3];
1640 AliESDtrack ctrack(*track); // Propagate on copy track
1641 if(!ctrack.PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
1643 Double_t pull[2] = {0, 0};
1644 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1645 for(Int_t i=0; i<2; i++){
1646 if(error[i]!=0) pull[i] = dz[i]/error[i]; // unitless
1650 if(charge<0) iPart = 0; // electron
1651 if(charge>0) iPart = 1; // positron
1653 printf("this is not an electron! Check HFEpid method\n");
1657 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1658 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]) {
1659 fHistHPDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1660 fHistHPDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1661 fHistHPDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1662 fHistHPDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1668 fHistDataHfePid[iPart]->Fill(esdpt);