]> 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
Changes for #82873: Module debugging broken (Christian)
[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**************************************************************************/
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"
faee3b18 32#include "AliMCVertex.h"
33
34#include "AliKFParticle.h"
35#include "AliKFVertex.h"
36
37#include "AliESDVertex.h"
38
39#include "AliPID.h"
70da6c5a 40
41#include "AliHFEdca.h"
42
faee3b18 43
70da6c5a 44ClassImp(AliHFEdca)
45
46//________________________________________________________________________
47const Char_t* AliHFEdca::fgkParticles[12] = {
faee3b18 48 // particles name
49 "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
50 "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
51 "allNegative", "allPositive"
70da6c5a 52};
faee3b18 53
54const 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
70da6c5a 60//________________________________________________________________________
61const Int_t AliHFEdca::fgkColorPart[12] = {
faee3b18 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
70da6c5a 66};
67
70da6c5a 68//________________________________________________________________________
faee3b18 69const 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};
70da6c5a 77
78//________________________________________________________________________
79const Char_t *AliHFEdca::fgkDcaVar[2] = {
faee3b18 80 "DcaXY", "DcaZ"};
70da6c5a 81
82//________________________________________________________________________
83const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
faee3b18 84 ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
85
86//________________________________________________________________________
87const Char_t *AliHFEdca::fgkVertexVar[3] = {
88 "VertexX", "VertexY", "VertexZ"};
89
90//________________________________________________________________________
91const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
92 ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
93
94//________________________________________________________________________
95const Char_t *AliHFEdca::fgkResDcaVar[2] = {
96 "deltaDcaXY", "deltaDcaZ"};
70da6c5a 97
faee3b18 98//________________________________________________________________________
99const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
100 ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
70da6c5a 101
102//________________________________________________________________________
103const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
faee3b18 104 "pullDcaXY", "pullDcaZ"
70da6c5a 105};
106
107//________________________________________________________________________
108const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
faee3b18 109 ";residual dca_{xy}/(error dca_{xy});counts",
110 ";residual dca_{z}/(error dca_{z});counts"
70da6c5a 111};
112
113//________________________________________________________________________
faee3b18 114const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
115 ";dca_{xy}^{data}/error dca_{xy};counts",
116 ";dca_{z}^{data}/error dca_{z};counts"
117};
70da6c5a 118
faee3b18 119//________________________________________________________________________
120AliHFEdca::AliHFEdca():
121 TObject(),
122 fStat(NULL)
70da6c5a 123{
faee3b18 124 // default constructor
125
69ac0e6f 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
70da6c5a 142}
143
144//________________________________________________________________________
faee3b18 145AliHFEdca::AliHFEdca(const AliHFEdca &ref):
146 TObject(ref),
147 fStat(ref.fStat)
70da6c5a 148{
faee3b18 149 // copy constructor
70da6c5a 150
69ac0e6f 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
70da6c5a 167}
168//_______________________________________________________________________________________________
faee3b18 169AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
70da6c5a 170{
faee3b18 171 //
172 // Assignment operator
173 //
174
faee3b18 175 if(this == &ref) return *this;
176 TObject::operator=(ref);
177 return *this;
69ac0e6f 178
70da6c5a 179}
180
181//________________________________________________________________________
182AliHFEdca::~AliHFEdca()
183{
faee3b18 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];
70da6c5a 196
faee3b18 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];
69ac0e6f 242
243
244 // Data
faee3b18 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 }
70da6c5a 257
faee3b18 258 if(fStat) delete fStat;
70da6c5a 259
6555e2ad 260 //Printf("analysis done\n");
69ac0e6f 261
70da6c5a 262}
263
264//________________________________________________________________________
c2690925 265void AliHFEdca::InitAnalysis()const{
faee3b18 266
6555e2ad 267 //Printf("initialize analysis\n");
70da6c5a 268
269}
270
271
272//________________________________________________________________________
273void AliHFEdca::PostAnalysis() const
274{
faee3b18 275 // do fit
276 // moved to dcaPostAnalysis.C
70da6c5a 277
278}
faee3b18 279
280
70da6c5a 281//________________________________________________________________________
282void AliHFEdca::CreateHistogramsResidual(TList *residualList){
faee3b18 283 // define histogram
284 // 1. residual
285
286 // for residuals
69ac0e6f 287 // fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
288 // fHistDcaZRes[kNParticles][kNPtBins]=0x0;
289 // fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
290 // fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
faee3b18 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++){
70da6c5a 302
303 TString histName((const char*)fgkParticles[j]);
faee3b18 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]);
70da6c5a 307
308 if(k==0){
309 fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
6555e2ad 310 fHistDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 311 if(j<(kNParticles-2)){
312 fHistEPDcaXYRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 313 fHistEPDcaXYRes[j][i]->SetLineColor((int)fgkColorPart[j]);}
70da6c5a 314 }
315 if(k==1){
316 fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 317 fHistDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 318 if(j<(kNParticles-2)){
319 fHistEPDcaZRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 320 fHistEPDcaZRes[j][i]->SetLineColor((int)fgkColorPart[j]); }
70da6c5a 321 }
faee3b18 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
70da6c5a 341
342}
343
344
345//________________________________________________________________________
346void AliHFEdca::CreateHistogramsPull(TList *pullList){
faee3b18 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 -----------------------------------------------------------------------
69ac0e6f 356 // fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
357 // fHistDcaZPull[kNParticles][kNPtBins]=0x0;
358 // fHistEPDcaXYPull[kNParticles-2][kNPtBins]=0x0;
359 // fHistEPDcaZPull[kNParticles-2][kNPtBins]=0x0;
faee3b18 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++){
70da6c5a 367
faee3b18 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]);
70da6c5a 372
373 if(k==0){
374 fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 375 fHistDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 376 if(j<(kNParticles-2)) {
377 fHistEPDcaXYPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 378 fHistEPDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
70da6c5a 379 }
380 if(k==1){
381 fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 382 fHistDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 383 if(j<(kNParticles-2)) {
384 fHistEPDcaZPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 385 fHistEPDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 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//________________________________________________________________________
408void 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
69ac0e6f 421 // fHistDcaXY[kNParticles][kNPtBins]=0x0;
422 // fHistDcaZ[kNParticles][kNPtBins]=0x0;
423 // fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
424 // fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
faee3b18 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);
69ac0e6f 445 fHistDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 446
447 if(j<(kNParticles-2)){
448 fHistEPDcaXY[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 449 fHistEPDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 450 }
451 if(k==1){
452 fHistDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 453 fHistDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 454 if(j<(kNParticles-2)){
455 fHistEPDcaZ[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 456 fHistEPDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);}
faee3b18 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//________________________________________________________________________
480void 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
69ac0e6f 493 // fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
494 // fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
faee3b18 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);
69ac0e6f 511 fHistKFDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 512 }
513 if(k==1){
514 fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 515 fHistKFDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 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//________________________________________________________________________
536void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
537 //
538 // define histograms: real Data
539 //
540
541 // for dca
69ac0e6f 542// fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
543// fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
544// fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
545// fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
faee3b18 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);
69ac0e6f 564 fHistDataDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 565
566 fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
69ac0e6f 567 fHistDataWoDcaXY[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 568 }
569 if(k==1){
570 fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 571 fHistDataDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 572
573 fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
69ac0e6f 574 fHistDataWoDcaZ[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 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//________________________________________________________________________
598void 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 -----------------------------------------------------------------------
69ac0e6f 607// fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
608// fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
faee3b18 609
69ac0e6f 610// fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
611// fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
faee3b18 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);
69ac0e6f 628 fHistDataDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 629
630 fHistDataWoDcaXYPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
69ac0e6f 631 fHistDataWoDcaXYPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 632 }
633 if(k==1){
634 fHistDataDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 635 fHistDataDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
faee3b18 636
637 fHistDataWoDcaZPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
69ac0e6f 638 fHistDataWoDcaZPull[j][i]->SetLineColor((int)fgkColorPart[j]);
70da6c5a 639 }
faee3b18 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//________________________________________________________________________
660void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
661 //
662 // define histograms: vertex
663 //
664 // for vertex
665
69ac0e6f 666// fHistMCvertex[kNVertexVar]=0x0;
667// fHistESDvertex[kNVertexVar]=0x0;
faee3b18 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//________________________________________________________________________
707void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
708 //
709 // define histograms: vertex
710 //
711 // for data vertex
712
69ac0e6f 713// fHistDatavertex[kNVertexVar]=0x0;
faee3b18 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//_______________________________________________________________________________________________
746void 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//_______________________________________________________________________________________________
780void 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//_______________________________________________________________________________________________
c2690925 809void AliHFEdca::FillHistogramsDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, AliMCEvent * const mcEvent)
faee3b18 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
c2690925 842 AliESDtrack ctrack(*track);
843 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
70da6c5a 844
69ac0e6f 845 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
846 if(!mctrack) return;
faee3b18 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//_______________________________________________________________________________________________
c2690925 940void AliHFEdca::FillHistogramsKfDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliMCEvent * const mcEvent)
6555e2ad 941 {
faee3b18 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];
c2690925 960 AliESDtrack ctrack(*track);
961 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 962
963 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
69ac0e6f 964 if(!mctrack) return;
faee3b18 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//_______________________________________________________________________________________________
c2690925 1060void AliHFEdca::FillHistogramsVtx(const AliESDEvent *const esdEvent, const AliMCEvent *const mcEvent)
faee3b18 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//_______________________________________________________________________________________________
c2690925 1085void AliHFEdca::FillHistogramsPid(const AliESDtrack * const track, const AliMCEvent * const mcEvent)
faee3b18 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())));
69ac0e6f 1095 if(!mctrack) return;
faee3b18 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]);
70da6c5a 1125}
1126
1127
faee3b18 1128////_______________________________________________________________________________________________
c2690925 1129void AliHFEdca::FillHistogramsDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
faee3b18 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
c2690925 1154 AliESDtrack ctrack(*track);
1155 if(!ctrack.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};
c2690925 1168 if(!ctrack.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//_______________________________________________________________________________________________
c2690925 1220void AliHFEdca::FillHistogramsDataVtx(const AliESDEvent * const esdEvent)
faee3b18 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////_______________________________________________________________________________________________
c2690925 1241void AliHFEdca::FillHistogramsDataPid(const 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//_________________________________________________________________________________________________
c2690925 1266void AliHFEdca::ApplyExtraCuts(const AliESDEvent * const esdEvent, Int_t nMinPrimVtxContributor)
faee3b18 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//_____________________________________________________
c2690925 1281Int_t AliHFEdca::GetCombinedPid(const AliESDtrack *const track)
faee3b18 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//_______________________________________________________________________________________________
c2690925 1508void AliHFEdca::FillHistogramsHfeDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliMCEvent * const mcEvent)
faee3b18 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];
c2690925 1539 AliESDtrack ctrack(*track);
1540 if(!ctrack.PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 1541
1542 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
69ac0e6f 1543 if(!mctrack) return;
faee3b18 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//_______________________________________________________________________________________________
c2690925 1616void AliHFEdca::FillHistogramsHfeDataDca(const AliESDEvent * const esdEvent, const AliESDtrack * const track, const AliESDVertex * const vtxESDSkip)
faee3b18 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
70da6c5a 1626// obtaining errors of dca ------------------------------------------------------------------
6555e2ad 1627 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex(); // UNUSED!
faee3b18 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.;
6555e2ad 1636
faee3b18 1637 Double_t dz[2]; // error of dca in cm
1638 Double_t covardz[3];
6555e2ad 1639
c2690925 1640 AliESDtrack ctrack(*track); // Propagate on copy track
1641 if(!ctrack.PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
faee3b18 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++){
6555e2ad 1646 if(error[i]!=0) pull[i] = dz[i]/error[i]; // unitless
faee3b18 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 }
6555e2ad 1665 else continue;
faee3b18 1666 }
1667
1668 fHistDataHfePid[iPart]->Fill(esdpt);
1669
70da6c5a 1670}
1671