]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEdca.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possibility to select on the basis of the presence of at least a fake track among...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEdca.cxx
CommitLineData
70da6c5a 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**************************************************************************/
27de2dfb 15
16/* $Id$ */
17
70da6c5a 18//
19// Class for impact parameter (DCA) of charged particles
20// Study resolution and pull: prepare for beauty study
21//
22// Authors:
23// Hongyan Yang <hongyan@physi.uni-heidelberg.de>
24// Carlo Bombonati <carlo.bombonati@cern.ch>
25//
26
27#include "TMath.h"
28#include "TH1F.h"
29#include "TList.h"
30#include <TParticle.h>
31#include "AliMCParticle.h"
32#include "AliESDtrack.h"
33#include "AliESDEvent.h"
34#include "AliMCEvent.h"
faee3b18 35#include "AliMCVertex.h"
36
37#include "AliKFParticle.h"
38#include "AliKFVertex.h"
39
40#include "AliESDVertex.h"
41
42#include "AliPID.h"
70da6c5a 43
44#include "AliHFEdca.h"
45
faee3b18 46
70da6c5a 47ClassImp(AliHFEdca)
48
49//________________________________________________________________________
50const Char_t* AliHFEdca::fgkParticles[12] = {
faee3b18 51 // particles name
52 "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
53 "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
54 "allNegative", "allPositive"
70da6c5a 55};
faee3b18 56
57const Int_t AliHFEdca::fgkPdgParticle[10] = {
58// 11, 13, -211, -233, -2122,
59// -11, -13, 211, 233, 2122};
60 kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
61 -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
62
70da6c5a 63//________________________________________________________________________
64const Int_t AliHFEdca::fgkColorPart[12] = {
faee3b18 65 // colors assigned to particles
66 kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
67 kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
68 kBlack, kGray+1
70da6c5a 69};
70
70da6c5a 71//________________________________________________________________________
faee3b18 72const Float_t AliHFEdca::fgkPtIntv[51] = {
73 // define pT bins
74 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55,
75 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50,
76 1.70, 1.90, 2.10, 2.30, 2.50, 2.70, 2.90, 3.10, 3.30, 3.50,
77 3.80, 4.10, 4.40, 4.70, 5.00, 5.30, 5.60, 5.90, 6.20, 6.50,
78 7.00, 7.50, 8.00, 9.00, 10.0, 11.0, 12.0, 13.0, 15.0, 18.0,
79 20.00};
70da6c5a 80
81//________________________________________________________________________
82const Char_t *AliHFEdca::fgkDcaVar[2] = {
faee3b18 83 "DcaXY", "DcaZ"};
70da6c5a 84
85//________________________________________________________________________
86const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
faee3b18 87 ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
88
89//________________________________________________________________________
90const Char_t *AliHFEdca::fgkVertexVar[3] = {
91 "VertexX", "VertexY", "VertexZ"};
92
93//________________________________________________________________________
94const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
95 ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
96
97//________________________________________________________________________
98const Char_t *AliHFEdca::fgkResDcaVar[2] = {
99 "deltaDcaXY", "deltaDcaZ"};
70da6c5a 100
faee3b18 101//________________________________________________________________________
102const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
103 ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
70da6c5a 104
105//________________________________________________________________________
106const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
faee3b18 107 "pullDcaXY", "pullDcaZ"
70da6c5a 108};
109
110//________________________________________________________________________
111const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
faee3b18 112 ";residual dca_{xy}/(error dca_{xy});counts",
113 ";residual dca_{z}/(error dca_{z});counts"
70da6c5a 114};
115
116//________________________________________________________________________
faee3b18 117const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
118 ";dca_{xy}^{data}/error dca_{xy};counts",
119 ";dca_{z}^{data}/error dca_{z};counts"
120};
70da6c5a 121
faee3b18 122//________________________________________________________________________
123AliHFEdca::AliHFEdca():
124 TObject(),
125 fStat(NULL)
70da6c5a 126{
faee3b18 127 // default constructor
128
69ac0e6f 129 for(Int_t j=0; j<kNParticles; j++){
130 fHistMcPid[j] = new TH1F();
131 fHistEsdPid[j] = new TH1F();
132 fHistDataEsdPid[j] = new TH1F();
133 }
134
135 for(Int_t i=0; i<3; i++){
136 fHistMCvertex[i] = new TH1F();
137 fHistESDvertex[i] = new TH1F();
138 fHistDatavertex[i] = new TH1F();
139 }
140
141 for(Int_t iEle=0; iEle<2; iEle++){
142 fHistDataHfePid[iEle] = new TH1F();
143 }
144
70da6c5a 145}
146
147//________________________________________________________________________
faee3b18 148AliHFEdca::AliHFEdca(const AliHFEdca &ref):
149 TObject(ref),
150 fStat(ref.fStat)
70da6c5a 151{
faee3b18 152 // copy constructor
70da6c5a 153
69ac0e6f 154 for(Int_t j=0; j<kNParticles; j++){
155 fHistMcPid[j] = ref.fHistMcPid[j];
156 fHistEsdPid[j] = ref.fHistEsdPid[j];
157 fHistDataEsdPid[j] = ref.fHistDataEsdPid[j];
158 }
159
160 for(Int_t i=0; i<3; i++){
161 fHistMCvertex[i] = ref.fHistMCvertex[i];
162 fHistESDvertex[i] = ref.fHistESDvertex[i];
163 fHistDatavertex[i] = ref.fHistDatavertex[i];
164 }
165
166 for(Int_t iEle=0; iEle<2; iEle++){
167 fHistDataHfePid[iEle] = ref.fHistDataHfePid[iEle];
168 }
169
70da6c5a 170}
171//_______________________________________________________________________________________________
faee3b18 172AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
70da6c5a 173{
faee3b18 174 //
175 // Assignment operator
176 //
177
faee3b18 178 if(this == &ref) return *this;
179 TObject::operator=(ref);
180 return *this;
69ac0e6f 181
70da6c5a 182}
183
184//________________________________________________________________________
185AliHFEdca::~AliHFEdca()
186{
faee3b18 187 // default destructor
188
189 for(Int_t j=0; j<kNParticles; j++){
190 for(Int_t i=0; i<kNPtBins; i++){
191 if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
192 if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
193
194 if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
195 if(fHistDcaZPull[j][i]) delete fHistDcaZPull[j][i];
196
197 if(fHistDcaXY[j][i]) delete fHistDcaXY[j][i];
198 if(fHistDcaZ[j][i]) delete fHistDcaZ[j][i];
70da6c5a 199
faee3b18 200 if(j<(kNParticles-2)){
201 if(fHistEPDcaXYRes[j][i]) delete fHistEPDcaXYRes[j][i];
202 if(fHistEPDcaZRes[j][i]) delete fHistEPDcaZRes[j][i];
203
204 if(fHistEPDcaXYPull[j][i]) delete fHistEPDcaXYPull[j][i];
205 if(fHistEPDcaZPull[j][i]) delete fHistEPDcaZPull[j][i];
206
207 if(fHistEPDcaXY[j][i]) delete fHistEPDcaXY[j][i];
208 if(fHistEPDcaZ[j][i]) delete fHistEPDcaZ[j][i];
209 }
210
211 if(fHistKFDcaXY[j][i]) delete fHistKFDcaXY[j][i];
212 if(fHistKFDcaZ[j][i]) delete fHistKFDcaZ[j][i];
213
214 if(fHistDataDcaXY[j][i]) delete fHistDataDcaXY[j][i];
215 if(fHistDataDcaZ[j][i]) delete fHistDataDcaZ[j][i];
216 if(fHistDataWoDcaXY[j][i]) delete fHistDataWoDcaXY[j][i];
217 if(fHistDataWoDcaZ[j][i]) delete fHistDataWoDcaZ[j][i];
218
219 if(fHistDataDcaXYPull[j][i]) delete fHistDataDcaXYPull[j][i];
220 if(fHistDataDcaZPull[j][i]) delete fHistDataDcaZPull[j][i];
221 if(fHistDataWoDcaXYPull[j][i]) delete fHistDataWoDcaXYPull[j][i];
222 if(fHistDataWoDcaZPull[j][i]) delete fHistDataWoDcaZPull[j][i];
223 }
224
225 if(fHistMcPid[j]) delete fHistMcPid[j];
226 if(fHistEsdPid[j]) delete fHistEsdPid[j];
227 if(fHistDataEsdPid[j]) delete fHistDataEsdPid[j];
228 }
229
230 for(Int_t i=0; i<3; i++){
231 if(fHistMCvertex[i]) delete fHistMCvertex[i];
232 if(fHistESDvertex[i]) delete fHistESDvertex[i];
233 if(fHistDatavertex[i]) delete fHistDatavertex[i];
234 }
235
236 // for the HFEpid
237 for(Int_t iEle=0; iEle<2; iEle++){
238 for(Int_t iPt=0; iPt<kNPtBins; iPt++){
239 if(fHistHPDcaXYRes[iEle][iPt]) delete fHistHPDcaXYRes[iEle][iPt];
240 if(fHistHPDcaZRes[iEle][iPt]) delete fHistHPDcaZRes[iEle][iPt];
241 if(fHistHPDcaXYPull[iEle][iPt]) delete fHistHPDcaXYPull[iEle][iPt];
242 if(fHistHPDcaZPull[iEle][iPt]) delete fHistHPDcaZPull[iEle][iPt];
243 if(fHistHPDcaXY[iEle][iPt]) delete fHistHPDcaXY[iEle][iPt];
244 if(fHistHPDcaZ[iEle][iPt]) delete fHistHPDcaZ[iEle][iPt];
69ac0e6f 245
246
247 // Data
faee3b18 248 if(fHistHPDataDcaXY[iEle][iPt]) delete fHistHPDataDcaXY[iEle][iPt];
249 if(fHistHPDataDcaZ[iEle][iPt]) delete fHistHPDataDcaZ[iEle][iPt];
250 if(fHistHPDataDcaXYPull[iEle][iPt]) delete fHistHPDataDcaXYPull[iEle][iPt];
251 if(fHistHPDataDcaZPull[iEle][iPt]) delete fHistHPDataDcaZPull[iEle][iPt];
252
253 }
254 for(Int_t i=0; i<2; i++)
255 if(fHistHfePid[iEle][i]) delete fHistHfePid[iEle][i];
256
257 if(fHistDataHfePid[iEle]) delete fHistDataHfePid[iEle];
258
259 }
70da6c5a 260
faee3b18 261 if(fStat) delete fStat;
70da6c5a 262
6555e2ad 263 //Printf("analysis done\n");
69ac0e6f 264
70da6c5a 265}
266
267//________________________________________________________________________
268void AliHFEdca::InitAnalysis(){
faee3b18 269
6555e2ad 270 //Printf("initialize analysis\n");
70da6c5a 271
272}
273
274
275//________________________________________________________________________
276void AliHFEdca::PostAnalysis() const
277{
faee3b18 278 // do fit
279 // moved to dcaPostAnalysis.C
70da6c5a 280
281}
faee3b18 282
283
70da6c5a 284//________________________________________________________________________
285void AliHFEdca::CreateHistogramsResidual(TList *residualList){
faee3b18 286 // define histogram
287 // 1. residual
288
289 // for residuals
69ac0e6f 290 // fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
291 // fHistDcaZRes[kNParticles][kNPtBins]=0x0;
292 // fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
293 // fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
faee3b18 294
295 const Int_t nBins = 1000;
296 const Float_t maxXYBin = 1000.;
297 const Float_t maxZBin = 1000.;
298
299
300 for(Int_t k=0; k<kNDcaVar; k++){
301 TString histTitle((const char*)fgkResDcaVarTitle[k]);
302
303 for(Int_t j=0; j<kNParticles; j++){
304 for(Int_t i=0; i<kNPtBins; i++){
70da6c5a 305
306 TString histName((const char*)fgkParticles[j]);
faee3b18 307 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
308 TString histEPName((const char*)fgkParticles[j]);
309 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
70da6c5a 310
311 if(k==0){
312 fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
6555e2ad 313 fHistDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 314 if(j<(kNParticles-2)){
315 fHistEPDcaXYRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 316 fHistEPDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);}
70da6c5a 317 }
318 if(k==1){
319 fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 320 fHistDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 321 if(j<(kNParticles-2)){
322 fHistEPDcaZRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 323 fHistEPDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]); }
70da6c5a 324 }
faee3b18 325 } // 50 pt bins
326 } //12 nparticles
327 } // 2 dca var
328
329 // TList *fResidualList = 0;
330 residualList->SetOwner();
331 residualList->SetName("residual");
332 for(Int_t iPart=0; iPart<kNParticles; iPart++){
333 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
334 residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
335 residualList->Add(fHistDcaZRes[iPart][iPtBin]);
336 if(iPart<(kNParticles-2)){
337 residualList->Add(fHistEPDcaXYRes[iPart][iPtBin]);
338 residualList->Add(fHistEPDcaZRes[iPart][iPtBin]);
339 }
340 } // loop over pt bins
341 } // loop over particles (pos, neg)
342
343
70da6c5a 344
345}
346
347
348//________________________________________________________________________
349void AliHFEdca::CreateHistogramsPull(TList *pullList){
faee3b18 350 // define histogram
351 // 2. pull
352
353 const Int_t nBins = 1000;
354 const Float_t maxXYBin = 20.;
355 const Float_t maxZBin = 20.;
356
357
358 // for pull -----------------------------------------------------------------------
69ac0e6f 359 // fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
360 // fHistDcaZPull[kNParticles][kNPtBins]=0x0;
361 // fHistEPDcaXYPull[kNParticles-2][kNPtBins]=0x0;
362 // fHistEPDcaZPull[kNParticles-2][kNPtBins]=0x0;
faee3b18 363
364
365 for(Int_t k=0; k<kNDcaVar; k++){
366 TString histTitle((const char*)fgkPullDcaVarTitle[k]);
367
368 for(Int_t j=0; j<kNParticles; j++){
369 for(Int_t i=0; i<kNPtBins; i++){
70da6c5a 370
faee3b18 371 TString histName((const char*)fgkParticles[j]);
372 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
373 TString histEPName((const char*)fgkParticles[j]);
374 histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
70da6c5a 375
376 if(k==0){
377 fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 378 fHistDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 379 if(j<(kNParticles-2)) {
380 fHistEPDcaXYPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 381 fHistEPDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
70da6c5a 382 }
383 if(k==1){
384 fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 385 fHistDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 386 if(j<(kNParticles-2)) {
387 fHistEPDcaZPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 388 fHistEPDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 389 }
390 } // 50 pt bins
391 } //6 nparticles
392 } // 2 dca var
393
394 // TList *fPullList = 0;
395 pullList->SetOwner();
396 pullList->SetName("pull");
397 for(Int_t iPart=0; iPart<kNParticles; iPart++){
398 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
399 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
400 pullList->Add(fHistDcaZPull[iPart][iPtBin]);
401 if(iPart<(kNParticles-2)){
402 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
403 pullList->Add(fHistDcaZPull[iPart][iPtBin]); }
404 } // loop over pt bins
405 } // loop over particles (pos, neg)
406
407}
408
409
410//________________________________________________________________________
411void AliHFEdca::CreateHistogramsDca(TList *dcaList){
412 //
413 // define histograms: MC dca
414 //
415
416 // statistics
417 fStat = 0x0;
418 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
419 fStat->SetMarkerStyle(20);
420 fStat->SetMarkerColor(3);
421 fStat->SetMarkerSize(1);
422
423 // for dca
69ac0e6f 424 // fHistDcaXY[kNParticles][kNPtBins]=0x0;
425 // fHistDcaZ[kNParticles][kNPtBins]=0x0;
426 // fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
427 // fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
faee3b18 428
429 const Int_t nBins = 1000;
430 const Float_t maxXYBin = 1000.;
431 const Float_t maxZBin = 1000.;
432
433
434 for(Int_t k=0; k<kNDcaVar; k++){
435 TString histTitle((const char*)fgkDcaVarTitle[k]);
436
437 for(Int_t j=0; j<kNParticles; j++){
438 for(Int_t i=0; i<kNPtBins; i++){
439
440 TString histName((const char*)fgkParticles[j]);
441 histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
442
443 TString histNameEP((const char*)fgkParticles[j]);
444 histNameEP += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
445
446 if(k==0){
447 fHistDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 448 fHistDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 449
450 if(j<(kNParticles-2)){
451 fHistEPDcaXY[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 452 fHistEPDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 453 }
454 if(k==1){
455 fHistDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 456 fHistDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 457 if(j<(kNParticles-2)){
458 fHistEPDcaZ[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 459 fHistEPDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 460 }
461 } // 50 pt bins
462 } //12 nparticles
463 } // 2 dca var
464
465 // TList *fDcaList = 0;
466 dcaList->SetOwner();
467 dcaList->SetName("mcDcaDistr");
468 for(Int_t iPart=0; iPart<kNParticles; iPart++){
469 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
470 dcaList->Add(fHistDcaXY[iPart][iPtBin]);
471 dcaList->Add(fHistDcaZ[iPart][iPtBin]);
472 if(iPart<(kNParticles-2)) {
473 dcaList->Add(fHistEPDcaXY[iPart][iPtBin]);
474 dcaList->Add(fHistEPDcaZ[iPart][iPtBin]); }
475 } // loop over pt bins
476 } // loop over particles (pos, neg)
477
478 dcaList->Add(fStat);
479
480}
481
482//________________________________________________________________________
483void AliHFEdca::CreateHistogramsKfDca(TList *kfDcaList){
484 //
485 // define histograms: MC dca
486 //
487
488 // statistics
489 fStat = 0x0;
490 fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
491 fStat->SetMarkerStyle(20);
492 fStat->SetMarkerColor(3);
493 fStat->SetMarkerSize(1);
494
495 // for kf dca
69ac0e6f 496 // fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
497 // fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
faee3b18 498
499 const Int_t nBins = 1000;
500 const Float_t maxXYBin = 1000.;
501 const Float_t maxZBin = 1000.;
502
503
504 for(Int_t k=0; k<kNDcaVar; k++){
505 TString histTitle((const char*)fgkDcaVarTitle[k]);
506
507 for(Int_t j=0; j<kNParticles; j++){
508 for(Int_t i=0; i<kNPtBins; i++){
509 TString histNameKF((const char*)fgkParticles[j]);
510 histNameKF += Form("_MCpid_KF%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
511
512 if(k==0){
513 fHistKFDcaXY[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 514 fHistKFDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 515 }
516 if(k==1){
517 fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 518 fHistKFDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 519 }
520 } // 50 pt bins
521 } //12 nparticles
522 } // 2 dca var
523
524 kfDcaList->SetOwner();
525 kfDcaList->SetName("mcKfDcaDistr");
526 for(Int_t iPart=0; iPart<kNParticles; iPart++){
527 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
528 kfDcaList->Add(fHistKFDcaXY[iPart][iPtBin]);
529 kfDcaList->Add(fHistKFDcaZ[iPart][iPtBin]);
530 } // loop over pt bins
531 } // loop over particles (pos, neg)
532
533 kfDcaList->Add(fStat);
534
535}
536
537
538//________________________________________________________________________
539void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
540 //
541 // define histograms: real Data
542 //
543
544 // for dca
69ac0e6f 545// fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
546// fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
547// fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
548// fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
faee3b18 549
550 const Int_t nBins = 1000;
551 const Float_t maxXYBin = 1000.;
552 const Float_t maxZBin = 1000.;
553
554 for(Int_t k=0; k<kNDcaVar; k++){
555 TString histTitle((const char*)fgkDcaVarTitle[k]);
556 for(Int_t j=0; j<kNParticles; j++){
557 for(Int_t i=0; i<kNPtBins; i++){
558
559 TString histName((const char*)fgkParticles[j]);
560 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
561
562 TString histNameWo((const char*)fgkParticles[j]);
563 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
564
565 if(k==0){
566 fHistDataDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 567 fHistDataDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 568
569 fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 570 fHistDataWoDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 571 }
572 if(k==1){
573 fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 574 fHistDataDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 575
576 fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 577 fHistDataWoDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 578 }
579 } // 50 pt bins
580 } //12 nparticles
581 } // 2 dca var
582
583 // TList *fDcaList = 0;
584 dataDcaList->SetOwner();
585 dataDcaList->SetName("dataDcaDistr");
586 for(Int_t iPart=0; iPart<kNParticles; iPart++){
587 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
588 dataDcaList->Add(fHistDataDcaXY[iPart][iPtBin]);
589 dataDcaList->Add(fHistDataDcaZ[iPart][iPtBin]);
590
591 dataDcaList->Add(fHistDataWoDcaXY[iPart][iPtBin]);
592 dataDcaList->Add(fHistDataWoDcaZ[iPart][iPtBin]);
593 } // loop over pt bins
594 } // loop over particles (pos, neg)
595
596
597}
598
599
600//________________________________________________________________________
601void AliHFEdca::CreateHistogramsDataPull(TList *dataPullList){
602 // define histogram
603 // 2. pull
604
605 const Int_t nBins = 1000;
606 const Float_t maxXYBin = 20.;
607 const Float_t maxZBin = 20.;
608
609 // for pull -----------------------------------------------------------------------
69ac0e6f 610// fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
611// fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
faee3b18 612
69ac0e6f 613// fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
614// fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
faee3b18 615
616
617 for(Int_t k=0; k<kNDcaVar; k++){
618 TString histTitle((const char*)fgkPullDataDcaVarTitle[k]);
619
620 for(Int_t j=0; j<kNParticles; j++){
621 for(Int_t i=0; i<kNPtBins; i++){
622
623 TString histName((const char*)fgkParticles[j]);
624 histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
625
626 TString histNameWo((const char*)fgkParticles[j]);
627 histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
628
629 if(k==0){
630 fHistDataDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 631 fHistDataDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 632
633 fHistDataWoDcaXYPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 634 fHistDataWoDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 635 }
636 if(k==1){
637 fHistDataDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 638 fHistDataDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 639
640 fHistDataWoDcaZPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 641 fHistDataWoDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
70da6c5a 642 }
faee3b18 643 } // 50 pt bins
644 } //6 nparticles
645 } // 2 dca var
646
647 // TList *fDataPullList = 0;
648 dataPullList->SetOwner();
649 dataPullList->SetName("dataPull");
650 for(Int_t iPart=0; iPart<kNParticles; iPart++){
651 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
652 dataPullList->Add(fHistDataDcaXYPull[iPart][iPtBin]);
653 dataPullList->Add(fHistDataDcaZPull[iPart][iPtBin]);
654
655 dataPullList->Add(fHistDataWoDcaXYPull[iPart][iPtBin]);
656 dataPullList->Add(fHistDataWoDcaZPull[iPart][iPtBin]);
657 } // loop over pt bins
658 } // loop over particles (pos, neg)
659
660}
661
662//________________________________________________________________________
663void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
664 //
665 // define histograms: vertex
666 //
667 // for vertex
668
69ac0e6f 669// fHistMCvertex[kNVertexVar]=0x0;
670// fHistESDvertex[kNVertexVar]=0x0;
faee3b18 671
672 const Int_t nBins = 1000;
673 const Float_t minXBin = -0.2e4;
674 const Float_t maxXBin = 0.2e4;
675 const Float_t minYBin = -0.5e4;
676 const Float_t maxYBin = 0.5e4;
677 const Float_t minZBin = -1.5e5;
678 const Float_t maxZBin = 1.5e5;
679
680 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
681 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
682
683 for(Int_t k=0; k<kNVertexVar; k++){
684 TString histTitle((const char*)fgkVertexVarTitle[k]);
685 TString histNameMC((const char*)fgkVertexVar[k]);
686 histNameMC += Form("_MC");
687 TString histNameESD((const char*)fgkVertexVar[k]);
688 histNameESD += Form("_ESD");
689
690 fHistMCvertex[k] = new TH1F((const char*)histNameMC, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
691 fHistMCvertex[k]->SetLineColor(k+2);
692
693 fHistESDvertex[k] = new TH1F((const char*)histNameESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
694 fHistESDvertex[k]->SetLineColor(k+2);
695 } // 3 vertex var
696
697 vertexList->SetOwner();
698 vertexList->SetName("vertexDistr");
699
700 for(Int_t k=0; k<kNVertexVar; k++){
701 vertexList->Add(fHistMCvertex[k]);
702 vertexList->Add(fHistESDvertex[k]);
703 }
704
705}
706
707
708
709//________________________________________________________________________
710void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
711 //
712 // define histograms: vertex
713 //
714 // for data vertex
715
69ac0e6f 716// fHistDatavertex[kNVertexVar]=0x0;
faee3b18 717
718 const Int_t nBins = 1000;
719 const Float_t minXBin = -0.2e4;
720 const Float_t maxXBin = 0.2e4;
721 const Float_t minYBin = -0.5e4;
722 const Float_t maxYBin = 0.5e4;
723 const Float_t minZBin = -1.5e5;
724 const Float_t maxZBin = 1.5e5;
725
726 const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
727 const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
728
729 for(Int_t k=0; k<kNVertexVar; k++){
730 TString histTitle((const char*)fgkVertexVarTitle[k]);
731 TString histNameDataESD((const char*)fgkVertexVar[k]);
732 histNameDataESD += Form("_data");
733
734 fHistDatavertex[k] = new TH1F((const char*)histNameDataESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
735 fHistDatavertex[k]->SetLineColor(k+2);
736 } // 3 vertex var
737
738 // TList *fVDaraVertexList = 0;
739 dataVertexList->SetOwner();
740 dataVertexList->SetName("dataVertexDistr");
741
742 for(Int_t k=0; k<kNVertexVar; k++){
743 dataVertexList->Add(fHistDatavertex[k]);
744 }
745
746}
747
748//_______________________________________________________________________________________________
749void AliHFEdca::CreateHistogramsPid(TList *mcPidList){
750 //
751 // define histograms which fills combined PID
752 //
753
754 const Char_t *mcOResd[2]={"mcPt", "esdPt"};
755
756 for(Int_t iPart=0; iPart<kNParticles; iPart++){
757 TString histTitleMc((const char*)fgkParticles[iPart]);
758 TString histTitleEsd((const char*)fgkParticles[iPart]);
759 histTitleMc += Form("_McPid_%s;p_{T} [GeV/c];counts", mcOResd[0]);
760 histTitleEsd += Form("_EsdPid_%s;p_{T} [GeV/c];counts", mcOResd[1]);
761
762 TString histNameMc((const char*)fgkParticles[iPart]);
763 TString histNameEsd((const char*)fgkParticles[iPart]);
764 histNameMc+=Form("_McPid_%s", mcOResd[0]);
765 histNameEsd+=Form("_EsdPid_%s", mcOResd[1]);
766
767 fHistMcPid[iPart] = new TH1F(histNameMc, histTitleMc, kNPtBins, fgkPtIntv);
768 fHistEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
769 }
770
771
772 mcPidList->SetOwner();
773 mcPidList->SetName("combinedPid");
774
775 for(Int_t iPart=0; iPart<kNParticles; iPart++){
776 mcPidList->Add(fHistMcPid[iPart]);
777 mcPidList->Add(fHistEsdPid[iPart]);
778 }
779}
780
781
782//_______________________________________________________________________________________________
783void AliHFEdca::CreateHistogramsDataPid(TList *pidList){
784 //
785 // define histograms which fills combined PID: data
786 //
787
788
789
790 for(Int_t iPart=0; iPart<kNParticles; iPart++){
791 TString histTitleEsd((const char*)fgkParticles[iPart]);
792 histTitleEsd+=Form("_DataEsdPid_esdPt;p_{T} [GeV/c];counts");
793 TString histNameEsd((const char*)fgkParticles[iPart]);
794 histNameEsd+=Form("_DataEsdPid");
795
796 fHistDataEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
797 }
798
799
800 pidList->SetOwner();
801 pidList->SetName("dataCombinedPid");
802
803 for(Int_t iPart=0; iPart<kNParticles; iPart++)
804 pidList->Add(fHistDataEsdPid[iPart]);
805
806
807}
808
809
810
811//_______________________________________________________________________________________________
812void AliHFEdca::FillHistogramsDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
813{
814 // the kDca plugin
815 // MC vertex
816 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
817 Double_t mcPrimV[3];
818 mcPrimV[0] = mcPrimVtx->GetX();
819 mcPrimV[1] = mcPrimVtx->GetY();
820 mcPrimV[2] = mcPrimVtx->GetZ();
821
822 Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
823
824// filling historgams track by track
825// obtaining reconstructed dca ------------------------------------------------------------------
826
827 Float_t esdpx = track->Px();
828 Float_t esdpy = track->Py();
829 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
830
831// obtaining errors of dca ------------------------------------------------------------------
832 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
833 Double_t primV[3];
834 primV[0] = primVtx->GetXv();
835 primV[1] = primVtx->GetYv();
836 primV[2] = primVtx->GetZv();
837
838 Float_t magneticField = 0; // initialized as 5kG
839 magneticField = esdEvent->GetMagneticField(); // in kG
840
841 Double_t beampiperadius=3.;
842 Double_t dz[2]; // error of dca in cm
843 Double_t covardz[3];
844
845 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
70da6c5a 846
69ac0e6f 847 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
848 if(!mctrack) return;
faee3b18 849 TParticle *part = mctrack->Particle();
850
851 Float_t vx = part->Vx(); // in cm
852 Float_t vy = part->Vy(); // in cm
853 Float_t vz = part->Vz(); // in cm
854
855 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
856
857 Float_t mcpx = part->Px();
858 Float_t mcpy = part->Py();
859 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
860
861 Int_t pdg = part->GetPdgCode();
862 Int_t esdPid = GetCombinedPid(track);
863
864 Int_t charge = 1;
865 if(pdg==kPDGelectron || pdg==kPDGmuon
866 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
867
868 // calculate mcDca ------------------------------------------------------------------
869 const Float_t conv[2] = {1.783/1.6, 2.99792458};
870 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
871
872 Float_t nx = esdpx/mcpt;
873 Float_t ny = esdpy/mcpt;
874
875 Float_t radius;
876 radius = TMath::Abs(radiusMc);
877
878 Double_t dxy = vxy - mcVtxXY; // in cm
879 Double_t dvx = vx - mcPrimV[0]; // in cm
880 Double_t dvy = vy - mcPrimV[1]; // in cm
881
882 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
883
884 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
885 Double_t residual[2] = {0, 0};
886 Double_t pull[2] = {0, 0};
887 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
888 for(Int_t i=0; i<2; i++){
889 residual[i] = dz[i] - mcDca[i]; // in centimeters
890 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
891 }
892
893
894 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
895 // identified ones
896 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
897 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
898 if(pdg==fgkPdgParticle[iPart]) {
899 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
900 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
901 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
902 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
903 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
904 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
905 } // mc pdg
906
907 if(esdPid==fgkPdgParticle[iPart]) {
908 fHistEPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
909 fHistEPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
910 fHistEPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
911 fHistEPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
912 fHistEPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
913 fHistEPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
914 } // esd pid
915
916 } // pt range
917
918 else
919 continue;
920 } // pt loop
921 } // particle id loop
922
923 // for charged particles: no pid
924 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
925 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
926 Int_t iPart = 10;
927 if(charge>0) iPart = 11;
928 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
929 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
930 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
931 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
932 fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
933 fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
934 }
935 else
936 continue;
937 } // pt
938
939}
940
941//_______________________________________________________________________________________________
6555e2ad 942void AliHFEdca::FillHistogramsKfDca(AliESDEvent * const esdEvent, AliESDtrack * const track, const AliMCEvent * const mcEvent)
943 {
faee3b18 944 // the kKfDca plugin
945
946// filling historgams track by track
947
948// obtaining reconstructed dca ------------------------------------------------------------------
949 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
950 Float_t magneticField = 0; // initialized as 5kG
951 magneticField = esdEvent->GetMagneticField(); // in kG
952
953 Float_t esdpx = track->Px();
954 Float_t esdpy = track->Py();
955 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
956
957 Int_t charge = track->Charge();
958
959 Double_t beampiperadius=3.;
960 Double_t dz[2]; // error of dca in cm
961 Double_t covardz[3];
962 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 963
964 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
69ac0e6f 965 if(!mctrack) return;
faee3b18 966 TParticle *part = mctrack->Particle();
967 Int_t pdg = part->GetPdgCode();
968
969 // calculate dca using AliKFParticle class------------------------------------------------------------------
970 Double_t kfdz[3] = {0, 0, 0};
971 Double_t kfdzwith[3] = {0, 0, 0};
972
973 Int_t trkID = track->GetID();
974
975 AliKFParticle::SetField(magneticField);
976 AliKFParticle kfParticle(*track, pdg);
977
978 // prepare kfprimary vertex
979 AliKFVertex kfESDprimary;
980 // Reconstruct Primary Vertex (with ESD tracks)
981 Int_t n=primVtx->GetNIndices();
982 if (n>0 && primVtx->GetStatus()){
983 kfESDprimary = AliKFVertex(*primVtx);
984
985 Double_t dcaXYWithTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
986 Double_t dcaWithTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
987 Double_t dcaZWithTrk = 0;
988 if(TMath::Abs(dcaWithTrk)>=TMath::Abs(dcaXYWithTrk))
989 dcaZWithTrk =TMath::Sqrt(dcaWithTrk*dcaWithTrk-dcaXYWithTrk*dcaXYWithTrk)*((dz[1]*-1<=0)?1:-1);
990 kfdzwith[0] = dcaXYWithTrk;
991 kfdzwith[1] = dcaZWithTrk;
992 kfdzwith[2] = dcaWithTrk; // with current track
993
994 Double_t dcaXYWoTrk = 0;
995 Double_t dcaZWoTrk = 0;
996 Double_t dcaWoTrk = 0;
997
998 UShort_t *priIndex = primVtx->GetIndices();
999
1000 for (Int_t i=0;i<n;i++){
1001
1002 Int_t idx = Int_t(priIndex[i]);
1003 if (idx == trkID){
1004 kfESDprimary -= kfParticle;
1005 dcaXYWoTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
1006 dcaWoTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
1007 if((dcaWoTrk-dcaXYWoTrk)>=0)
1008 dcaZWoTrk = TMath::Abs(dcaWoTrk*dcaWoTrk - dcaXYWoTrk*dcaXYWoTrk)*((dz[1]*-1<=0)?1:-1);
1009 } // remove current track from this calculation
1010 } // loop over all primary vertex contributors
1011
1012
1013 kfdz[0] = dcaXYWoTrk;
1014 kfdz[1] = dcaZWoTrk;
1015 kfdz[2] = dcaWoTrk;
1016
1017 } // only if n contributor > 0 and primVtx constructed
1018
1019 fStat->Fill(0);
1020
1021 if(dz[0]!=0 && dz[0]*kfdzwith[0]>0 && TMath::Abs(kfdzwith[0]/dz[0])>0.9999 && TMath::Abs(kfdzwith[0]/dz[0])<1.0001)fStat->Fill(1);; // same
1022 if(dz[0]!=0 && dz[0]*kfdzwith[0]<0 && TMath::Abs(kfdzwith[0]/dz[0])>0.9999 && TMath::Abs(kfdzwith[0]/dz[0])<1.0001) fStat->Fill(2); // swapped sign
1023 if(kfdzwith[0]==0 && dz[0]!=0) fStat->Fill(3); // 0 from KF particle (with current track)
1024
1025 if(dz[0]!=0 && dz[0]*kfdz[0]>0 && TMath::Abs(kfdz[0]/dz[0])>0.8 && TMath::Abs(kfdz[0]/dz[0])<1.2) fStat->Fill(-1);; // same
1026 if(dz[0]!=0 && dz[0]*kfdz[0]<0 && TMath::Abs(kfdz[0]/dz[0])>0.8 && TMath::Abs(kfdz[0]/dz[0]) <1.2) fStat->Fill(-2); // swapped sign
1027 if(kfdz[0]==0 && dz[0]!=0) fStat->Fill(-3); // 0 from KF particle (without current track)
1028
1029 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1030 // identified ones
1031 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1032 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1033 if(pdg==fgkPdgParticle[iPart]) {
1034 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1.0e4);
1035 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1.0e4);
1036 } // mc pdg
1037 } // pt range
1038
1039 else
1040 continue;
1041 } // pt loop
1042 } // particle id loop
1043
1044 // for charged particles: no pid
1045 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1046 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1047 Int_t iPart = 10;
1048 if(charge>0) iPart = 11;
1049 fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1050 fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1051
1052 }
1053 else
1054 continue;
1055 } // pt
1056
1057} // KF dca
1058
1059
1060//_______________________________________________________________________________________________
1061void AliHFEdca::FillHistogramsVtx(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
1062{
1063
1064 // MC vertex
1065 AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1066 Double_t mcPrimV[3];
1067 mcPrimV[0] = mcPrimVtx->GetX();
1068 mcPrimV[1] = mcPrimVtx->GetY();
1069 mcPrimV[2] = mcPrimVtx->GetZ();
1070
1071// obtaining errors of dca ------------------------------------------------------------------
1072 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1073 Double_t primV[3];
1074 primV[0] = primVtx->GetXv();
1075 primV[1] = primVtx->GetYv();
1076 primV[2] = primVtx->GetZv();
1077
1078 for(Int_t i=0; i<kNVertexVar; i++){
1079 fHistMCvertex[i]->Fill(mcPrimV[i]*1.0e4);
1080 fHistESDvertex[i]->Fill(primV[i]*1.0e4);
1081 }
1082
1083}
1084
1085//_______________________________________________________________________________________________
6555e2ad 1086void AliHFEdca::FillHistogramsPid(AliESDtrack * const track, const AliMCEvent * const mcEvent)
faee3b18 1087{
1088
1089
1090// filling historgams track by track
1091 Float_t esdpx = track->Px();
1092 Float_t esdpy = track->Py();
1093 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1094
1095 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
69ac0e6f 1096 if(!mctrack) return;
faee3b18 1097 TParticle *part = mctrack->Particle();
1098
1099 Float_t mcpx = part->Px();
1100 Float_t mcpy = part->Py();
1101 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1102
1103 Int_t pdg = part->GetPdgCode();
1104 Int_t esdPid = GetCombinedPid(track);
1105
1106
1107 Double_t ptMom[2] = {mcpt, esdpt};
1108 // for combined PID
1109 for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1110 if(pdg==fgkPdgParticle[iPart]) // pid all by MC
1111 fHistMcPid[iPart]->Fill(ptMom[0]);
1112
1113 if(esdPid==fgkPdgParticle[iPart]) // pid all by combined pid
1114 fHistEsdPid[iPart]->Fill(ptMom[1]);
1115 } // loop over particles
1116
1117 // for charged
1118 if(pdg==kPDGelectron || pdg==kPDGmuon || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton)
1119 fHistMcPid[10]->Fill(ptMom[0]);
1120 if(pdg==-kPDGelectron || pdg==-kPDGmuon || pdg==kPDGpion || pdg==kPDGkaon || pdg==kPDGproton)
1121 fHistMcPid[11]->Fill(ptMom[0]);
1122 if(esdPid==kPDGelectron || esdPid==kPDGmuon || esdPid==-kPDGpion || esdPid==-kPDGkaon || esdPid==-kPDGproton)
1123 fHistEsdPid[10]->Fill(ptMom[1]);
1124 if(esdPid==-kPDGelectron || esdPid==-kPDGmuon || esdPid==kPDGpion || esdPid==kPDGkaon || esdPid==kPDGproton)
1125 fHistEsdPid[11]->Fill(ptMom[1]);
70da6c5a 1126}
1127
1128
faee3b18 1129////_______________________________________________________________________________________________
1130void AliHFEdca::FillHistogramsDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1131{
1132// filling historgams track by track
1133// obtaining reconstructed dca --------------------------------------------------------------
1134
1135 Float_t esdpx = track->Px();
1136 Float_t esdpy = track->Py();
1137 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1138 Int_t charge = track->Charge();
1139
1140// obtaining errors of dca ------------------------------------------------------------------
1141 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1142 Double_t primV[3];
1143 primV[0] = primVtx->GetXv();
1144 primV[1] = primVtx->GetYv();
1145 primV[2] = primVtx->GetZv();
1146
1147
1148 Float_t magneticField = 0; // initialized as 5kG
1149 magneticField = esdEvent->GetMagneticField(); // in kG
1150
1151 Double_t beampiperadius=3.;
1152 Double_t dz[2]; // error of dca in cm
1153 Double_t covardz[3];
1154
1155 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
6555e2ad 1156
faee3b18 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(!track->PropagateToDCA(vtxESDSkip, magneticField, beampiperadius, dzwo, covardzwo)) return; // protection
faee3b18 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
70da6c5a 1219//_______________________________________________________________________________________________
faee3b18 1220void AliHFEdca::FillHistogramsDataVtx(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////_______________________________________________________________________________________________
1241void AliHFEdca::FillHistogramsDataPid(AliESDtrack * const track)
70da6c5a 1242{
1243// filling historgams track by track
faee3b18 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
70da6c5a 1257
faee3b18 1258 else {
1259 if(charge<0) fHistDataEsdPid[10]->Fill(esdpt);
1260 if(charge>0) fHistDataEsdPid[11]->Fill(esdpt);
1261 }
1262 }
1263}
1264
1265//_________________________________________________________________________________________________
1266void AliHFEdca::ApplyExtraCuts(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//_____________________________________________________
1281Int_t AliHFEdca::GetCombinedPid(AliESDtrack *const track)
1282{
1283
69ac0e6f 1284 // combined detector pid
faee3b18 1285 Double_t prob[AliPID::kSPECIES];
1286 track->GetESDpid(prob);
69ac0e6f 1287
1288 // setting priors!
1289 Double_t priors[AliPID::kSPECIESN];
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;
faee3b18 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//________________________________________________________________________
1344void 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
69ac0e6f 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;
faee3b18 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);
69ac0e6f 1384 fHistHPDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1385 fHistHPDcaXYRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
69ac0e6f 1386 fHistHPDcaXYRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1387 fHistHPDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
69ac0e6f 1388 fHistHPDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1389 }
1390
1391 if(k==1){
1392 fHistHPDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
69ac0e6f 1393 fHistHPDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1394 fHistHPDcaZRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
69ac0e6f 1395 fHistHPDcaZRes[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1396 fHistHPDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxZBinPull, 1+maxZBinPull);
69ac0e6f 1397 fHistHPDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1398 }
1399 } // 50 pt bins
1400 } //2 nparticles
1401 } // 2 dca var
1402
69ac0e6f 1403 // fHistHfePid[2][2] = 0x0; //! HFE pid
faee3b18 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//________________________________________________________________________
1433void 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
69ac0e6f 1447// fHistHPDataDcaXY[2][kNPtBins]=0x0;
1448// fHistHPDataDcaZ[2][kNPtBins]=0x0;
1449// fHistHPDataDcaXYPull[2][kNPtBins]=0x0;
1450// fHistHPDataDcaZPull[2][kNPtBins]=0x0;
faee3b18 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);
69ac0e6f 1464 fHistHPDataDcaXY[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1465 fHistHPDataDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
69ac0e6f 1466 fHistHPDataDcaXYPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1467 }
1468
1469 if(k==1){
1470 fHistHPDataDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
69ac0e6f 1471 fHistHPDataDcaZ[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1472 fHistHPDataDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsDca, 1-maxZBinPull, 1+maxZBinPull);
69ac0e6f 1473 fHistHPDataDcaZPull[iPart][i]->SetLineColor((int)fgkColorPart[iPart*5]);
faee3b18 1474 }
1475
1476 } // 50 pt bins
1477 } // 2 particle type
1478 } // 2 dca var
1479
6555e2ad 1480 //fHistDataHfePid[2] = 0x0; //! HFE pid
faee3b18 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//_______________________________________________________________________________________________
1508void AliHFEdca::FillHistogramsHfeDca(AliESDEvent * const esdEvent, AliESDtrack * const track, 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
70da6c5a 1521// obtaining reconstructed dca ------------------------------------------------------------------
faee3b18 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 if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 1540
1541 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
69ac0e6f 1542 if(!mctrack) return;
faee3b18 1543 TParticle *part = mctrack->Particle();
1544
1545 Float_t vx = part->Vx(); // in cm
1546 Float_t vy = part->Vy(); // in cm
1547 Float_t vz = part->Vz(); // in cm
1548
1549 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1550
1551 Float_t mcpx = part->Px();
1552 Float_t mcpy = part->Py();
1553 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1554
1555 Int_t pdg = part->GetPdgCode();
1556
1557 Int_t charge = 1;
1558 if(pdg==kPDGelectron || pdg==kPDGmuon
1559 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1560
1561 // calculate mcDca ------------------------------------------------------------------
1562 const Float_t conv[2] = {1.783/1.6, 2.99792458};
1563 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
1564
1565 Float_t nx = esdpx/mcpt;
1566 Float_t ny = esdpy/mcpt;
1567
1568 Float_t radius;
1569 radius = TMath::Abs(radiusMc);
1570
1571 Double_t dxy = vxy - mcVtxXY; // in cm
1572 Double_t dvx = vx - mcPrimV[0]; // in cm
1573 Double_t dvy = vy - mcPrimV[1]; // in cm
1574
1575 Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ; // in meters
1576
1577 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
1578 Double_t residual[2] = {0, 0};
1579 Double_t pull[2] = {0, 0};
1580 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1581 for(Int_t i=0; i<2; i++){
1582 residual[i] = dz[i] - mcDca[i]; // in centimeters
1583 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
1584 }
1585
1586 Int_t iPart = -1;
1587 if(track->Charge()<0) iPart = 0; // electron
1588 if(track->Charge()>0) iPart = 1; // positron
1589 if(track->Charge()==0) {
1590 printf("this is not an electron! Check HFEpid method");
1591 return;
1592 }
1593 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1594 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1595 fHistHPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
1596 fHistHPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
1597 fHistHPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1598 fHistHPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1599 fHistHPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
1600 fHistHPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
1601
1602 } // pt range
1603
1604 else
1605 continue;
1606 } // pt loop
1607
1608 fHistHfePid[iPart][0]->Fill(esdpt);
1609 fHistHfePid[iPart][1]->Fill(mcpt);
1610
1611}
1612
1613
1614//_______________________________________________________________________________________________
6555e2ad 1615void AliHFEdca::FillHistogramsHfeDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
faee3b18 1616{
1617// filling historgams track by track
1618// obtaining reconstructed dca --------------------------------------------------------------
1619
1620 Float_t esdpx = track->Px();
1621 Float_t esdpy = track->Py();
1622 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1623 Int_t charge = track->Charge();
1624
70da6c5a 1625// obtaining errors of dca ------------------------------------------------------------------
6555e2ad 1626 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!
faee3b18 1627 Double_t primV[3];
1628 primV[0] = primVtx->GetXv();
1629 primV[1] = primVtx->GetYv();
1630 primV[2] = primVtx->GetZv();
1631
1632 Float_t magneticField = 0; // initialized as 5kG
1633 magneticField = esdEvent->GetMagneticField(); // in kG
1634 Double_t beampiperadius=3.;
6555e2ad 1635
faee3b18 1636 Double_t dz[2]; // error of dca in cm
1637 Double_t covardz[3];
6555e2ad 1638
1639 if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 1640
1641 Double_t pull[2] = {0, 0};
1642 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1643 for(Int_t i=0; i<2; i++){
6555e2ad 1644 if(error[i]!=0) pull[i] = dz[i]/error[i]; // unitless
faee3b18 1645 }
1646
1647 Int_t iPart = -1;
1648 if(charge<0) iPart = 0; // electron
1649 if(charge>0) iPart = 1; // positron
1650 if(charge==0) {
1651 printf("this is not an electron! Check HFEpid method\n");
1652 return;
1653 }
1654
1655 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1656 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]) {
1657 fHistHPDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1658 fHistHPDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1659 fHistHPDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1660 fHistHPDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1661
1662 }
6555e2ad 1663 else continue;
faee3b18 1664 }
1665
1666 fHistDataHfePid[iPart]->Fill(esdpt);
1667
70da6c5a 1668}
1669