- Reshuffling of the particle codes in AliPID. Now the light nuclei are between the
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEdca.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
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 **************************************************************************/
15 //
16 // Class for impact parameter (DCA) of charged particles
17 // Study resolution and pull: prepare for beauty study
18 //
19 // Authors:
20 //   Hongyan Yang <hongyan@physi.uni-heidelberg.de>
21 //   Carlo Bombonati <carlo.bombonati@cern.ch>
22 //
23
24 #include "TMath.h"
25 #include "TH1F.h"
26 #include "TList.h"
27 #include <TParticle.h>
28 #include "AliMCParticle.h"
29 #include "AliESDtrack.h"
30 #include "AliESDEvent.h"
31 #include "AliMCEvent.h"
32 #include "AliMCVertex.h"
33
34 #include "AliKFParticle.h"
35 #include "AliKFVertex.h"
36
37 #include "AliESDVertex.h"
38
39 #include "AliPID.h"
40
41 #include "AliHFEdca.h"
42
43
44 ClassImp(AliHFEdca)
45
46 //________________________________________________________________________
47 const Char_t* AliHFEdca::fgkParticles[12] = {
48  // particles name
49  "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus", 
50  "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
51  "allNegative", "allPositive"
52 };
53
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};
59
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,
65  kBlack, kGray+1
66 };
67
68 //________________________________________________________________________
69 const Float_t AliHFEdca::fgkPtIntv[51] = {
70  // define pT bins
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, 
76  20.00};
77
78 //________________________________________________________________________
79 const Char_t *AliHFEdca::fgkDcaVar[2] = {
80  "DcaXY",  "DcaZ"};
81
82 //________________________________________________________________________
83 const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
84  ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
85
86 //________________________________________________________________________
87 const Char_t *AliHFEdca::fgkVertexVar[3] = {
88  "VertexX", "VertexY", "VertexZ"};
89
90 //________________________________________________________________________
91 const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
92  ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
93
94 //________________________________________________________________________
95 const Char_t *AliHFEdca::fgkResDcaVar[2] = {
96  "deltaDcaXY",  "deltaDcaZ"};
97
98 //________________________________________________________________________
99 const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
100  ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
101
102 //________________________________________________________________________
103 const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
104  "pullDcaXY", "pullDcaZ"
105 };
106
107 //________________________________________________________________________
108 const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
109  ";residual dca_{xy}/(error dca_{xy});counts", 
110  ";residual dca_{z}/(error dca_{z});counts"
111 };
112
113 //________________________________________________________________________
114 const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
115  ";dca_{xy}^{data}/error dca_{xy};counts", 
116  ";dca_{z}^{data}/error dca_{z};counts"
117 };
118
119 //________________________________________________________________________
120 AliHFEdca::AliHFEdca():
121   TObject(),
122   fStat(NULL)
123 {
124  // default constructor
125
126  for(Int_t j=0; j<kNParticles; j++){
127    fHistMcPid[j] = new TH1F();
128    fHistEsdPid[j] = new TH1F();
129    fHistDataEsdPid[j] = new TH1F();
130  }
131
132  for(Int_t i=0; i<3; i++){
133    fHistMCvertex[i] = new TH1F();
134    fHistESDvertex[i] = new TH1F();
135    fHistDatavertex[i] = new TH1F();   
136  }
137  
138  for(Int_t iEle=0; iEle<2; iEle++){
139    fHistDataHfePid[iEle] = new TH1F();
140  }
141
142 }
143
144 //________________________________________________________________________
145 AliHFEdca::AliHFEdca(const AliHFEdca &ref):
146  TObject(ref),
147  fStat(ref.fStat)
148 {
149  // copy constructor
150
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];
155  }
156
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];
161  }
162  
163  for(Int_t iEle=0; iEle<2; iEle++){
164    fHistDataHfePid[iEle] = ref.fHistDataHfePid[iEle];
165  }
166
167 }
168 //_______________________________________________________________________________________________
169 AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
170 {
171  //
172  // Assignment operator
173  //
174
175  if(this == &ref) return *this;
176  TObject::operator=(ref);
177  return *this;
178
179 }
180
181 //________________________________________________________________________
182 AliHFEdca::~AliHFEdca()
183 {
184  // default destructor
185
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];
190
191      if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
192      if(fHistDcaZPull[j][i]) delete fHistDcaZPull[j][i];
193
194      if(fHistDcaXY[j][i]) delete fHistDcaXY[j][i];
195      if(fHistDcaZ[j][i]) delete fHistDcaZ[j][i];
196
197      if(j<(kNParticles-2)){
198         if(fHistEPDcaXYRes[j][i]) delete fHistEPDcaXYRes[j][i];
199         if(fHistEPDcaZRes[j][i]) delete fHistEPDcaZRes[j][i];
200         
201         if(fHistEPDcaXYPull[j][i]) delete fHistEPDcaXYPull[j][i];
202         if(fHistEPDcaZPull[j][i]) delete fHistEPDcaZPull[j][i];
203         
204         if(fHistEPDcaXY[j][i]) delete fHistEPDcaXY[j][i];
205         if(fHistEPDcaZ[j][i]) delete fHistEPDcaZ[j][i];
206      }
207
208      if(fHistKFDcaXY[j][i]) delete fHistKFDcaXY[j][i];
209      if(fHistKFDcaZ[j][i]) delete fHistKFDcaZ[j][i];
210
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];
215
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];
220    }
221
222    if(fHistMcPid[j]) delete fHistMcPid[j];
223    if(fHistEsdPid[j]) delete fHistEsdPid[j];
224    if(fHistDataEsdPid[j]) delete fHistDataEsdPid[j];
225  }
226
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];    
231  }
232
233  // for the HFEpid
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];      
242      
243      
244      // Data
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];    
249
250    }
251    for(Int_t i=0; i<2; i++)
252      if(fHistHfePid[iEle][i]) delete fHistHfePid[iEle][i];
253
254    if(fHistDataHfePid[iEle]) delete fHistDataHfePid[iEle];
255
256  }
257
258  if(fStat) delete fStat;
259
260  //Printf("analysis done\n");
261
262 }
263
264 //________________________________________________________________________
265 void AliHFEdca::InitAnalysis()const{
266
267  //Printf("initialize analysis\n");
268
269 }
270
271
272 //________________________________________________________________________
273 void AliHFEdca::PostAnalysis() const
274 {
275  // do fit
276  // moved to dcaPostAnalysis.C
277
278 }
279
280
281 //________________________________________________________________________
282 void AliHFEdca::CreateHistogramsResidual(TList *residualList){
283  // define histogram
284  // 1. residual
285
286  // for residuals
287   // fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
288   // fHistDcaZRes[kNParticles][kNPtBins]=0x0;
289   // fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
290  // fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
291
292  const Int_t nBins = 1000;
293  const Float_t maxXYBin = 1000.;
294  const Float_t maxZBin = 1000.;
295
296
297  for(Int_t k=0; k<kNDcaVar; k++){
298    TString histTitle((const char*)fgkResDcaVarTitle[k]);
299
300    for(Int_t j=0; j<kNParticles; j++){
301      for(Int_t i=0; i<kNPtBins; i++){
302         
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]);
307         
308         if(k==0){
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]);}
314         }           
315         if(k==1){
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]); }
321         }   
322      } // 50 pt bins
323    } //12 nparticles
324  } // 2 dca var
325
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]);  
336      }
337    } // loop over pt bins
338  }  // loop over particles (pos, neg)
339
340
341
342 }
343
344
345 //________________________________________________________________________
346 void AliHFEdca::CreateHistogramsPull(TList *pullList){
347  // define histogram
348  // 2. pull
349
350  const Int_t nBins = 1000;
351  const Float_t maxXYBin = 20.;
352  const Float_t maxZBin = 20.;
353
354
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;
360
361
362  for(Int_t k=0; k<kNDcaVar; k++){
363    TString histTitle((const char*)fgkPullDcaVarTitle[k]);
364
365    for(Int_t j=0; j<kNParticles; j++){
366      for(Int_t i=0; i<kNPtBins; i++){
367         
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]);
372         
373         if(k==0){
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]);}
379         }           
380         if(k==1){
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]);}
386         }   
387      } // 50 pt bins
388    } //6 nparticles
389  } // 2 dca var
390
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)
403
404 }
405
406
407 //________________________________________________________________________
408 void AliHFEdca::CreateHistogramsDca(TList *dcaList){
409  // 
410  // define histograms: MC dca
411  //
412
413  // statistics
414  fStat = 0x0;
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); 
419
420  // for dca
421  // fHistDcaXY[kNParticles][kNPtBins]=0x0;
422  // fHistDcaZ[kNParticles][kNPtBins]=0x0;
423  // fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
424  // fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
425
426  const Int_t nBins = 1000;
427  const Float_t maxXYBin = 1000.;
428  const Float_t maxZBin = 1000.;
429
430
431  for(Int_t k=0; k<kNDcaVar; k++){
432    TString histTitle((const char*)fgkDcaVarTitle[k]);
433
434    for(Int_t j=0; j<kNParticles; j++){
435      for(Int_t i=0; i<kNPtBins; i++){
436         
437         TString histName((const char*)fgkParticles[j]); 
438         histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
439         
440         TString histNameEP((const char*)fgkParticles[j]);       
441         histNameEP += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
442
443         if(k==0){
444           fHistDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
445           fHistDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
446           
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]);}
450         }           
451         if(k==1){
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]);}
457         }   
458      } // 50 pt bins
459    } //12 nparticles
460  } // 2 dca var
461
462    //  TList *fDcaList = 0;
463  dcaList->SetOwner();
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)
474
475  dcaList->Add(fStat);
476
477 }
478
479 //________________________________________________________________________
480 void AliHFEdca::CreateHistogramsKfDca(TList *kfDcaList){
481  // 
482  // define histograms: MC dca
483  //
484
485  // statistics
486  fStat = 0x0;
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); 
491
492  // for kf dca
493  // fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
494  // fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
495
496  const Int_t nBins = 1000;
497  const Float_t maxXYBin = 1000.;
498  const Float_t maxZBin = 1000.;
499
500
501  for(Int_t k=0; k<kNDcaVar; k++){
502    TString histTitle((const char*)fgkDcaVarTitle[k]);
503
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]);
508         
509         if(k==0){
510           fHistKFDcaXY[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
511           fHistKFDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
512         }           
513         if(k==1){
514           fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
515           fHistKFDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
516         }   
517      } // 50 pt bins
518    } //12 nparticles
519  } // 2 dca var
520
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)
529
530  kfDcaList->Add(fStat);
531
532 }
533
534
535 //________________________________________________________________________
536 void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
537  //
538  // define histograms: real Data
539  //
540
541  // for dca
542 //  fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
543 //  fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
544 //  fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
545 //  fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
546
547  const Int_t nBins = 1000;
548  const Float_t maxXYBin = 1000.;
549  const Float_t maxZBin = 1000.;
550
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++){
555         
556         TString histName((const char*)fgkParticles[j]); 
557         histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
558
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]);
561         
562         if(k==0){
563           fHistDataDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
564           fHistDataDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
565
566           fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
567           fHistDataWoDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
568         }           
569         if(k==1){
570           fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
571           fHistDataDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
572
573           fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
574           fHistDataWoDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
575         }   
576      } // 50 pt bins
577    } //12 nparticles
578  } // 2 dca var
579
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]);  
587
588      dataDcaList->Add(fHistDataWoDcaXY[iPart][iPtBin]);  
589      dataDcaList->Add(fHistDataWoDcaZ[iPart][iPtBin]);  
590    } // loop over pt bins
591  }  // loop over particles (pos, neg)
592
593
594 }
595
596
597 //________________________________________________________________________
598 void AliHFEdca::CreateHistogramsDataPull(TList *dataPullList){
599  // define histogram
600  // 2. pull
601
602  const Int_t nBins = 1000;
603  const Float_t maxXYBin = 20.;
604  const Float_t maxZBin = 20.;
605
606  // for pull -----------------------------------------------------------------------
607 //  fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
608 //  fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
609
610 //  fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
611 //  fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
612
613
614  for(Int_t k=0; k<kNDcaVar; k++){
615    TString histTitle((const char*)fgkPullDataDcaVarTitle[k]);
616
617    for(Int_t j=0; j<kNParticles; j++){
618      for(Int_t i=0; i<kNPtBins; i++){
619         
620         TString histName((const char*)fgkParticles[j]); 
621         histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
622
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]);
625         
626         if(k==0){
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]);
629
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]);
632         }           
633         if(k==1){
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]);
636
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]);
639         }   
640      } // 50 pt bins
641    } //6 nparticles
642  } // 2 dca var
643
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]);  
651
652      dataPullList->Add(fHistDataWoDcaXYPull[iPart][iPtBin]);  
653      dataPullList->Add(fHistDataWoDcaZPull[iPart][iPtBin]);  
654    } // loop over pt bins
655  }  // loop over particles (pos, neg)
656
657 }
658
659 //________________________________________________________________________
660 void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
661  //
662  // define histograms: vertex
663  //
664  // for  vertex
665
666 //  fHistMCvertex[kNVertexVar]=0x0;
667 //  fHistESDvertex[kNVertexVar]=0x0;
668
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;
676
677  const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
678  const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
679
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");
686
687    fHistMCvertex[k] = new TH1F((const char*)histNameMC, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
688    fHistMCvertex[k]->SetLineColor(k+2);
689
690    fHistESDvertex[k] = new TH1F((const char*)histNameESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
691    fHistESDvertex[k]->SetLineColor(k+2);
692  } // 3 vertex var
693
694  vertexList->SetOwner();
695  vertexList->SetName("vertexDistr");
696
697  for(Int_t k=0; k<kNVertexVar; k++){
698    vertexList->Add(fHistMCvertex[k]);  
699    vertexList->Add(fHistESDvertex[k]);  
700  }
701
702 }
703
704
705
706 //________________________________________________________________________
707 void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
708  //
709  // define histograms: vertex
710  //
711  // for data vertex
712
713 //  fHistDatavertex[kNVertexVar]=0x0;
714
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;
722
723  const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
724  const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
725
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");
730
731    fHistDatavertex[k] = new TH1F((const char*)histNameDataESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
732    fHistDatavertex[k]->SetLineColor(k+2);
733  } // 3 vertex var
734
735    //  TList *fVDaraVertexList = 0;
736  dataVertexList->SetOwner();
737  dataVertexList->SetName("dataVertexDistr");
738
739  for(Int_t k=0; k<kNVertexVar; k++){
740    dataVertexList->Add(fHistDatavertex[k]);  
741  }
742
743 }
744
745 //_______________________________________________________________________________________________
746 void AliHFEdca::CreateHistogramsPid(TList *mcPidList){
747  //
748  // define histograms which fills combined PID
749  //
750
751  const Char_t *mcOResd[2]={"mcPt", "esdPt"};
752
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]);
758
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]);     
763
764    fHistMcPid[iPart] = new TH1F(histNameMc, histTitleMc, kNPtBins, fgkPtIntv);
765    fHistEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
766  }
767
768
769  mcPidList->SetOwner();
770  mcPidList->SetName("combinedPid");
771
772  for(Int_t iPart=0; iPart<kNParticles; iPart++){
773    mcPidList->Add(fHistMcPid[iPart]);
774    mcPidList->Add(fHistEsdPid[iPart]);
775  }
776 }
777
778
779 //_______________________________________________________________________________________________
780 void AliHFEdca::CreateHistogramsDataPid(TList *pidList){
781  //
782  // define histograms which fills combined PID: data
783  //
784
785
786
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");
792
793      fHistDataEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
794  }
795
796
797  pidList->SetOwner();
798  pidList->SetName("dataCombinedPid");
799
800  for(Int_t iPart=0; iPart<kNParticles; iPart++)
801    pidList->Add(fHistDataEsdPid[iPart]);
802
803
804 }
805
806
807
808 //_______________________________________________________________________________________________
809 void AliHFEdca::FillHistogramsDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, AliMCEvent * const mcEvent)
810 {
811  // the kDca plugin
812  // MC vertex
813  AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();      
814  Double_t mcPrimV[3];
815  mcPrimV[0] = mcPrimVtx->GetX();
816  mcPrimV[1] = mcPrimVtx->GetY();
817  mcPrimV[2] = mcPrimVtx->GetZ();
818
819  Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
820
821 // filling historgams track by track
822 // obtaining reconstructed dca ------------------------------------------------------------------
823
824  Float_t esdpx = track->Px();
825  Float_t esdpy = track->Py();
826  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);  
827
828 // obtaining errors of dca ------------------------------------------------------------------
829  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
830  Double_t primV[3];
831  primV[0] = primVtx->GetXv();
832  primV[1] = primVtx->GetYv();
833  primV[2] = primVtx->GetZv();
834
835  Float_t magneticField = 0;  // initialized as 5kG
836  magneticField = esdEvent->GetMagneticField();  // in kG
837
838  Double_t beampiperadius=3.;
839  Double_t dz[2];   // error of dca in cm
840  Double_t covardz[3];
841
842  AliESDtrack ctrack(*track);
843  if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return;  // protection
844
845  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))); 
846  if(!mctrack) return;
847  TParticle *part = mctrack->Particle();
848
849  Float_t vx = part->Vx();  // in cm
850  Float_t vy = part->Vy();  // in cm
851  Float_t vz = part->Vz();   // in cm
852
853  Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
854
855  Float_t mcpx = part->Px();
856  Float_t mcpy = part->Py();
857  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
858
859  Int_t pdg = part->GetPdgCode();
860  Int_t esdPid = GetCombinedPid(track);
861
862  Int_t charge = 1;
863  if(pdg==kPDGelectron || pdg==kPDGmuon 
864     || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;  
865
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
869
870  Float_t nx = esdpx/mcpt;
871  Float_t ny = esdpy/mcpt;
872
873  Float_t radius;
874  radius = TMath::Abs(radiusMc);
875
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
879
880  Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
881
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
889  }
890
891
892  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){    
893    // identified ones
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);  
903         }  // mc pdg      
904         
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);  
912         }  // esd pid
913         
914      } // pt range
915
916      else
917         continue;
918    }  // pt loop
919  } // particle id loop
920
921  // for charged particles: no pid
922  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
923    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
924      Int_t iPart = 10;
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);      
932    }
933    else
934      continue;
935  } // pt  
936
937 }
938
939 //_______________________________________________________________________________________________
940 void AliHFEdca::FillHistogramsKfDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliMCEvent * const mcEvent)
941  {
942  // the kKfDca plugin
943
944 // filling historgams track by track
945
946 // obtaining reconstructed dca ------------------------------------------------------------------  
947  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
948    Float_t magneticField = 0;  // initialized as 5kG
949  magneticField = esdEvent->GetMagneticField();  // in kG
950
951  Float_t esdpx = track->Px();
952  Float_t esdpy = track->Py();
953  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);  
954
955  Int_t charge = track->Charge();
956
957  Double_t beampiperadius=3.;  
958  Double_t dz[2];   // error of dca in cm
959  Double_t covardz[3];
960  AliESDtrack ctrack(*track);
961  if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection 
962
963  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));  
964  if(!mctrack) return;
965  TParticle *part = mctrack->Particle();    
966  Int_t pdg = part->GetPdgCode();
967
968  // calculate dca using AliKFParticle class------------------------------------------------------------------
969  Double_t kfdz[3] = {0, 0, 0};  
970  Double_t kfdzwith[3] = {0, 0, 0};  
971
972  Int_t trkID = track->GetID();
973
974  AliKFParticle::SetField(magneticField);
975  AliKFParticle kfParticle(*track, pdg);
976
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);
983
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
992
993    Double_t dcaXYWoTrk = 0;
994    Double_t dcaZWoTrk = 0;
995    Double_t dcaWoTrk = 0;
996
997    UShort_t *priIndex = primVtx->GetIndices();
998
999    for (Int_t i=0;i<n;i++){
1000
1001      Int_t idx = Int_t(priIndex[i]);
1002      if (idx == trkID){
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
1010
1011
1012    kfdz[0] = dcaXYWoTrk;  
1013    kfdz[1] = dcaZWoTrk;
1014    kfdz[2] = dcaWoTrk;  
1015
1016  }  // only if n contributor > 0 and primVtx constructed
1017
1018  fStat->Fill(0);
1019
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)
1023
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)
1027
1028  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){    
1029    // identified ones
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);  
1035         }  // mc pdg      
1036      } // pt range
1037
1038      else
1039         continue;
1040    }  // pt loop
1041  } // particle id loop
1042
1043  // for charged particles: no pid
1044  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1045    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1046      Int_t iPart = 10;
1047      if(charge>0) iPart = 11;
1048      fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1049      fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1050
1051    }
1052    else
1053      continue;
1054  } // pt  
1055
1056 }  // KF dca
1057
1058
1059 //_______________________________________________________________________________________________
1060 void AliHFEdca::FillHistogramsVtx(const AliESDEvent *const esdEvent, const AliMCEvent *const mcEvent)
1061 {
1062
1063  // MC vertex
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();
1069
1070 // obtaining errors of dca ------------------------------------------------------------------
1071  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
1072  Double_t primV[3];
1073  primV[0] = primVtx->GetXv();
1074  primV[1] = primVtx->GetYv();
1075  primV[2] = primVtx->GetZv();
1076
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); 
1080  }
1081
1082 }
1083
1084 //_______________________________________________________________________________________________
1085 void AliHFEdca::FillHistogramsPid(const AliESDtrack * const track, const AliMCEvent * const mcEvent)
1086 {
1087
1088
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);    
1093
1094  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));  
1095  if(!mctrack) return;
1096  TParticle *part = mctrack->Particle();
1097
1098  Float_t mcpx = part->Px();
1099  Float_t mcpy = part->Py();
1100  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1101
1102  Int_t pdg = part->GetPdgCode();
1103  Int_t esdPid = GetCombinedPid(track);
1104
1105
1106  Double_t ptMom[2] = {mcpt, esdpt};
1107  // for combined PID
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]);
1111
1112    if(esdPid==fgkPdgParticle[iPart])             // pid all by combined pid
1113      fHistEsdPid[iPart]->Fill(ptMom[1]);    
1114  } // loop over particles
1115
1116  // for charged
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]);  
1125 }
1126
1127
1128 ////_______________________________________________________________________________________________
1129 void AliHFEdca::FillHistogramsDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
1130 {
1131 // filling historgams track by track
1132 // obtaining reconstructed dca --------------------------------------------------------------
1133
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();
1138
1139 // obtaining errors of dca ------------------------------------------------------------------
1140  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
1141  Double_t primV[3];
1142  primV[0] = primVtx->GetXv();
1143  primV[1] = primVtx->GetYv();
1144  primV[2] = primVtx->GetZv();
1145
1146
1147  Float_t magneticField = 0;  // initialized as 5kG
1148  magneticField = esdEvent->GetMagneticField();  // in kG
1149
1150  Double_t beampiperadius=3.;
1151  Double_t dz[2];   // error of dca in cm
1152  Double_t covardz[3];
1153
1154  AliESDtrack ctrack(*track);
1155  if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return;  // protection
1156
1157
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                                                
1162  }
1163
1164  // get dca when current track is not included
1165
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
1169
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                                                
1173  }
1174
1175  // do pid 
1176  Int_t esdPid = GetCombinedPid(track);
1177
1178  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){    
1179    // identified ones
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]);    
1191      }
1192      else
1193         continue;
1194    }
1195  } 
1196
1197  // for charged particles
1198  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1199    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1200      Int_t iPart = 10;
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]);       
1211
1212    }
1213    else
1214      continue;
1215  } 
1216
1217 }
1218
1219 //_______________________________________________________________________________________________
1220 void AliHFEdca::FillHistogramsDataVtx(const AliESDEvent * const esdEvent)
1221 {
1222
1223
1224 // obtaining errors of dca ------------------------------------------------------------------
1225  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
1226  Double_t primV[3];
1227  primV[0] = primVtx->GetXv();
1228  primV[1] = primVtx->GetYv();
1229  primV[2] = primVtx->GetZv();
1230
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);
1236
1237 }
1238
1239
1240 ////_______________________________________________________________________________________________
1241 void AliHFEdca::FillHistogramsDataPid(const AliESDtrack * const track)
1242 {
1243 // filling historgams track by track
1244 // obtaining reconstructed dca --------------------------------------------------------------
1245
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();
1250
1251  Int_t esdPid = GetCombinedPid(track);
1252
1253  for(Int_t iPart=0; iPart<kNParticles; iPart++){
1254    if(iPart<(kNParticles-2)){
1255      if(esdPid==fgkPdgParticle[iPart])  fHistDataEsdPid[iPart]->Fill(esdpt);
1256    }     // for identified
1257
1258    else {
1259      if(charge<0) fHistDataEsdPid[10]->Fill(esdpt);
1260      if(charge>0) fHistDataEsdPid[11]->Fill(esdpt);
1261    }
1262  }
1263 }  
1264
1265 //_________________________________________________________________________________________________
1266 void AliHFEdca::ApplyExtraCuts(const AliESDEvent * const esdEvent, Int_t nMinPrimVtxContributor)
1267
1268
1269  //
1270  // only one extra cut, number of contributors to each primary vertex
1271  //
1272
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
1277
1278 }
1279
1280 //_____________________________________________________
1281 Int_t AliHFEdca::GetCombinedPid(const AliESDtrack *const track) 
1282 {
1283
1284  // combined detector pid             
1285  Double_t prob[AliPID::kSPECIES];
1286  track->GetESDpid(prob);
1287
1288  // setting priors!
1289  Double_t priors[AliPID::kSPECIES];
1290  priors[0] = 0.01;
1291  priors[1] = 0.01;
1292  priors[2] = 0.85;
1293  priors[3] = 0.10;
1294  priors[4] = 0.05;
1295
1296  Int_t charge = track->Charge();
1297  Int_t esdPid = -1; 
1298
1299  AliPID pid;
1300  pid.SetPriors(priors);
1301  pid.SetProbabilities(prob);
1302
1303  // identify particle as the most probable 
1304
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; 
1310
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; 
1316
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; 
1322
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; 
1328
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; 
1334
1335
1336  return charge*esdPid;
1337
1338 }
1339
1340
1341 // for the HFE pid
1342
1343 //________________________________________________________________________
1344 void AliHFEdca::CreateHistogramsHfeDca(TList *hfeDcaList){
1345  //
1346  // define histograms: hfe pid electrons in MC
1347  //
1348
1349  const Int_t nBinsDca = 1000;
1350  const Float_t maxXYBinDca = 1000.;
1351  const Float_t maxZBinDca = 1000.;
1352
1353  const Int_t nBinsPull = 1000;
1354  const Float_t maxXYBinPull = 20.;
1355  const Float_t maxZBinPull = 20.;
1356
1357  const Char_t *mcOResd[2]={"mcPt", "esdPt"};
1358
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;
1365
1366
1367  for(Int_t k=0; k<kNDcaVar; k++){
1368
1369    TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1370    TString histTitleRes((const char*)fgkPullDcaVarTitle[k]);
1371    TString histTitlePull((const char*)fgkResDcaVarTitle[k]);
1372
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]);
1381         
1382         if(k==0){
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]);
1389         }
1390         
1391         if(k==1){
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]);
1398         }
1399      } // 50 pt bins
1400    } //2 nparticles
1401  } // 2 dca var
1402
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);
1411    }
1412  }  
1413
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]);
1424    }
1425    for(Int_t id=0; id<2; id++)
1426      hfeDcaList->Add(fHistHfePid[id][iPart]);
1427  }
1428
1429 }
1430
1431
1432 //________________________________________________________________________
1433 void AliHFEdca::CreateHistogramsHfeDataDca(TList *hfeDataDcaList){
1434  //
1435  // define histograms: hfe pid electrons in data
1436  //
1437
1438  const Int_t nBinsDca = 1000;
1439  const Float_t maxXYBinDca = 1000.;
1440  const Float_t maxZBinDca = 1000.;
1441
1442  const Int_t nBinsPull = 1000;
1443  const Float_t maxXYBinPull = 20.;
1444  const Float_t maxZBinPull = 20.;
1445
1446
1447 //  fHistHPDataDcaXY[2][kNPtBins]=0x0;
1448 //  fHistHPDataDcaZ[2][kNPtBins]=0x0;
1449 //  fHistHPDataDcaXYPull[2][kNPtBins]=0x0;
1450 //  fHistHPDataDcaZPull[2][kNPtBins]=0x0;
1451
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]);
1461         
1462         if(k==0){
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]);
1467         }
1468         
1469         if(k==1){
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]);
1474         }
1475         
1476      } // 50 pt bins
1477    } // 2 particle type
1478  } // 2 dca var
1479
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);
1487  }
1488
1489
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]); 
1498
1499      hfeDataDcaList->Add(fHistDataHfePid[iPart]);
1500    }
1501  }
1502
1503
1504 }
1505
1506
1507 //_______________________________________________________________________________________________
1508 void AliHFEdca::FillHistogramsHfeDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliMCEvent * const mcEvent)
1509 {
1510  // the kHFEpid plugin
1511
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();
1517
1518  Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
1519
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);  
1525
1526 // obtaining errors of dca ------------------------------------------------------------------
1527  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();      
1528  Double_t primV[3];
1529  primV[0] = primVtx->GetXv();
1530  primV[1] = primVtx->GetYv();
1531  primV[2] = primVtx->GetZv();
1532
1533  Float_t magneticField = 0;  // initialized as 5kG
1534  magneticField = esdEvent->GetMagneticField();  // in kG
1535
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
1541
1542  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));  
1543  if(!mctrack) return;
1544  TParticle *part = mctrack->Particle();
1545
1546  Float_t vx = part->Vx();  // in cm
1547  Float_t vy = part->Vy();  // in cm
1548  Float_t vz = part->Vz();   // in cm
1549
1550  Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1551
1552  Float_t mcpx = part->Px();
1553  Float_t mcpy = part->Py();
1554  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1555
1556  Int_t pdg = part->GetPdgCode();
1557
1558  Int_t charge = 1;
1559  if(pdg==kPDGelectron || pdg==kPDGmuon 
1560     || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1561
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
1565
1566  Float_t nx = esdpx/mcpt;
1567  Float_t ny = esdpy/mcpt;
1568
1569  Float_t radius;
1570  radius = TMath::Abs(radiusMc);
1571
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
1575
1576  Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
1577
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
1585  }
1586
1587  Int_t iPart = -1;
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");
1592    return;
1593  }
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);  
1602
1603    } // pt range
1604
1605    else
1606      continue;
1607  }  // pt loop
1608
1609  fHistHfePid[iPart][0]->Fill(esdpt);
1610  fHistHfePid[iPart][1]->Fill(mcpt);
1611
1612 }
1613
1614
1615 //_______________________________________________________________________________________________
1616 void AliHFEdca::FillHistogramsHfeDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
1617 {
1618 // filling historgams track by track
1619 // obtaining reconstructed dca --------------------------------------------------------------
1620
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();
1625
1626 // obtaining errors of dca ------------------------------------------------------------------
1627  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!     
1628  Double_t primV[3];
1629  primV[0] = primVtx->GetXv();
1630  primV[1] = primVtx->GetYv();
1631  primV[2] = primVtx->GetZv();
1632
1633  Float_t magneticField = 0;  // initialized as 5kG
1634  magneticField = esdEvent->GetMagneticField();  // in kG
1635  Double_t beampiperadius=3.; 
1636
1637  Double_t dz[2];   // error of dca in cm
1638  Double_t covardz[3];
1639
1640  AliESDtrack ctrack(*track); // Propagate on copy track
1641  if(!ctrack.PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
1642
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    
1647  }
1648
1649  Int_t iPart = -1;
1650  if(charge<0) iPart = 0;  // electron
1651  if(charge>0) iPart = 1;  // positron
1652  if(charge==0) {
1653    printf("this is not an electron! Check HFEpid method\n");
1654    return;
1655  }
1656
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]);               
1663
1664    }
1665    else continue;
1666  }
1667
1668  fHistDataHfePid[iPart]->Fill(esdpt);
1669
1670 }
1671