PWG3/hfe/AliHFEdca.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Class for impact parameter (DCA) of charged particles
  17 // Study resolution and pull: prepare for beauty study
  18 //
  19 // Authors:
  20 //   Hongyan Yang <hongyan@physi.uni-heidelberg.de>
  21 //   Carlo Bombonati <carlo.bombonati@cern.ch>
  22 //
  23
  24 #include "TMath.h"
  25 #include "TH1F.h"
  26 #include "TList.h"
  27 #include <TParticle.h>
  28 #include "AliMCParticle.h"
  29 #include "AliESDtrack.h"
  30 #include "AliESDEvent.h"
  31 #include "AliMCEvent.h"
  32 #include "AliMCVertex.h"
  33
  34 #include "AliKFParticle.h"
  35 #include "AliKFVertex.h"
  36
  37 #include "AliESDVertex.h"
  38
  39 #include "AliPID.h"
  40
  41 #include "AliHFEdca.h"
  42
  43
  44 ClassImp(AliHFEdca)
  45
  46 //________________________________________________________________________
  47 const Char_t* AliHFEdca::fgkParticles[12] = {
  48  // particles name
  49  "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
  50  "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
  51  "allNegative", "allPositive"
  52 };
  53
  54 const Int_t AliHFEdca::fgkPdgParticle[10] = {
  55 //   11, 13, -211, -233, -2122,
  56 //   -11, -13, 211, 233, 2122};
  57  kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
  58  -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
  59
  60 //________________________________________________________________________
  61 const Int_t AliHFEdca::fgkColorPart[12] = {
  62  // colors assigned to particles
  63  kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
  64  kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
  65  kBlack, kGray+1
  66 };
  67
  68 //________________________________________________________________________
  69 const Float_t AliHFEdca::fgkPtIntv[51] = {
  70  // define pT bins
  71  0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55,
  72  0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50,
  73  1.70, 1.90, 2.10, 2.30, 2.50, 2.70, 2.90, 3.10, 3.30, 3.50,
  74  3.80, 4.10, 4.40, 4.70, 5.00, 5.30, 5.60, 5.90, 6.20, 6.50,
  75  7.00, 7.50, 8.00, 9.00, 10.0, 11.0, 12.0, 13.0, 15.0, 18.0,
  76  20.00};
  77
  78 //________________________________________________________________________
  79 const Char_t *AliHFEdca::fgkDcaVar[2] = {
  80  "DcaXY",  "DcaZ"};
  81
  82 //________________________________________________________________________
  83 const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
  84  ";dca_{xy} [#mum];counts", ";dca_{z} [#mum];counts"};
  85
  86 //________________________________________________________________________
  87 const Char_t *AliHFEdca::fgkVertexVar[3] = {
  88  "VertexX", "VertexY", "VertexZ"};
  89
  90 //________________________________________________________________________
  91 const Char_t *AliHFEdca::fgkVertexVarTitle[3] ={
  92  ";vertex_{x} [#mum];counts", ";vertex_{y} [#mum];counts", ";vertex_{z} [#mum];counts"};
  93
  94 //________________________________________________________________________
  95 const Char_t *AliHFEdca::fgkResDcaVar[2] = {
  96  "deltaDcaXY",  "deltaDcaZ"};
  97
  98 //________________________________________________________________________
  99 const Char_t *AliHFEdca::fgkResDcaVarTitle[2] ={
 100  ";residual #Delta(d_{xy}) [#mum];counts", ";residual #Delta(d_{z}) [#mum];counts"};
 101
 102 //________________________________________________________________________
 103 const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
 104  "pullDcaXY", "pullDcaZ"
 105 };
 106
 107 //________________________________________________________________________
 108 const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
 109  ";residual dca_{xy}/(error dca_{xy});counts",
 110  ";residual dca_{z}/(error dca_{z});counts"
 111 };
 112
 113 //________________________________________________________________________
 114 const Char_t *AliHFEdca::fgkPullDataDcaVarTitle[2] = {
 115  ";dca_{xy}^{data}/error dca_{xy};counts",
 116  ";dca_{z}^{data}/error dca_{z};counts"
 117 };
 118
 119 //________________________________________________________________________
 120 AliHFEdca::AliHFEdca():
 121   TObject(),
 122   fStat(NULL)
 123 {
 124  // default constructor
 125
 126 }
 127
 128 //________________________________________________________________________
 129 AliHFEdca::AliHFEdca(const AliHFEdca &ref):
 130  TObject(ref),
 131  fStat(ref.fStat)
 132 {
 133  // copy constructor
 134
 135 }
 136 //_______________________________________________________________________________________________
 137 AliHFEdca&AliHFEdca::operator=(const AliHFEdca &ref)
 138 {
 139  //
 140  // Assignment operator
 141  //
 142
 143
 144  if(this == &ref) return *this;
 145  TObject::operator=(ref);
 146  return *this;
 147 }
 148
 149 //________________________________________________________________________
 150 AliHFEdca::~AliHFEdca()
 151 {
 152  // default destructor
 153
 154  for(Int_t j=0; j<kNParticles; j++){
 155    for(Int_t i=0; i<kNPtBins; i++){
 156      if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
 157      if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
 158
 159      if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
 160      if(fHistDcaZPull[j][i]) delete fHistDcaZPull[j][i];
 161
 162      if(fHistDcaXY[j][i]) delete fHistDcaXY[j][i];
 163      if(fHistDcaZ[j][i]) delete fHistDcaZ[j][i];
 164
 165      if(j<(kNParticles-2)){
 166         if(fHistEPDcaXYRes[j][i]) delete fHistEPDcaXYRes[j][i];
 167         if(fHistEPDcaZRes[j][i]) delete fHistEPDcaZRes[j][i];
 168
 169         if(fHistEPDcaXYPull[j][i]) delete fHistEPDcaXYPull[j][i];
 170         if(fHistEPDcaZPull[j][i]) delete fHistEPDcaZPull[j][i];
 171
 172         if(fHistEPDcaXY[j][i]) delete fHistEPDcaXY[j][i];
 173         if(fHistEPDcaZ[j][i]) delete fHistEPDcaZ[j][i];
 174      }
 175
 176      if(fHistKFDcaXY[j][i]) delete fHistKFDcaXY[j][i];
 177      if(fHistKFDcaZ[j][i]) delete fHistKFDcaZ[j][i];
 178
 179      if(fHistDataDcaXY[j][i]) delete fHistDataDcaXY[j][i];
 180      if(fHistDataDcaZ[j][i]) delete fHistDataDcaZ[j][i];
 181      if(fHistDataWoDcaXY[j][i]) delete fHistDataWoDcaXY[j][i];
 182      if(fHistDataWoDcaZ[j][i]) delete fHistDataWoDcaZ[j][i];
 183
 184      if(fHistDataDcaXYPull[j][i]) delete fHistDataDcaXYPull[j][i];
 185      if(fHistDataDcaZPull[j][i]) delete fHistDataDcaZPull[j][i];
 186      if(fHistDataWoDcaXYPull[j][i]) delete fHistDataWoDcaXYPull[j][i];
 187      if(fHistDataWoDcaZPull[j][i]) delete fHistDataWoDcaZPull[j][i];
 188    }
 189
 190    if(fHistMcPid[j]) delete fHistMcPid[j];
 191    if(fHistEsdPid[j]) delete fHistEsdPid[j];
 192    if(fHistDataEsdPid[j]) delete fHistDataEsdPid[j];
 193  }
 194
 195  for(Int_t i=0; i<3; i++){
 196    if(fHistMCvertex[i]) delete fHistMCvertex[i];
 197    if(fHistESDvertex[i]) delete fHistESDvertex[i];
 198    if(fHistDatavertex[i]) delete fHistDatavertex[i];
 199  }
 200
 201  // for the HFEpid
 202  for(Int_t iEle=0; iEle<2; iEle++){
 203    for(Int_t iPt=0; iPt<kNPtBins; iPt++){
 204      if(fHistHPDcaXYRes[iEle][iPt]) delete fHistHPDcaXYRes[iEle][iPt];
 205      if(fHistHPDcaZRes[iEle][iPt]) delete fHistHPDcaZRes[iEle][iPt];
 206      if(fHistHPDcaXYPull[iEle][iPt]) delete fHistHPDcaXYPull[iEle][iPt];
 207      if(fHistHPDcaZPull[iEle][iPt]) delete fHistHPDcaZPull[iEle][iPt];
 208      if(fHistHPDcaXY[iEle][iPt]) delete fHistHPDcaXY[iEle][iPt];
 209      if(fHistHPDcaZ[iEle][iPt]) delete fHistHPDcaZ[iEle][iPt];
 210
 211
 212  // Data
 213      if(fHistHPDataDcaXY[iEle][iPt]) delete fHistHPDataDcaXY[iEle][iPt];
 214      if(fHistHPDataDcaZ[iEle][iPt]) delete fHistHPDataDcaZ[iEle][iPt];
 215      if(fHistHPDataDcaXYPull[iEle][iPt]) delete fHistHPDataDcaXYPull[iEle][iPt];
 216      if(fHistHPDataDcaZPull[iEle][iPt]) delete fHistHPDataDcaZPull[iEle][iPt];
 217
 218    }
 219    for(Int_t i=0; i<2; i++)
 220      if(fHistHfePid[iEle][i]) delete fHistHfePid[iEle][i];
 221
 222    if(fHistDataHfePid[iEle]) delete fHistDataHfePid[iEle];
 223
 224  }
 225
 226  if(fStat) delete fStat;
 227
 228  Printf("analysis done\n");
 229 }
 230
 231 //________________________________________________________________________
 232 void AliHFEdca::InitAnalysis(){
 233
 234  Printf("initialize analysis\n");
 235
 236 }
 237
 238
 239 //________________________________________________________________________
 240 void AliHFEdca::PostAnalysis() const
 241 {
 242  // do fit
 243  // moved to dcaPostAnalysis.C
 244
 245 }
 246
 247
 248 //________________________________________________________________________
 249 void AliHFEdca::CreateHistogramsResidual(TList *residualList){
 250  // define histogram
 251  // 1. residual
 252
 253  // for residuals
 254  fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
 255  fHistDcaZRes[kNParticles][kNPtBins]=0x0;
 256
 257  fHistEPDcaXYRes[kNParticles-2][kNPtBins]=0x0;
 258  fHistEPDcaZRes[kNParticles-2][kNPtBins]=0x0;
 259
 260  const Int_t nBins = 1000;
 261  const Float_t maxXYBin = 1000.;
 262  const Float_t maxZBin = 1000.;
 263
 264
 265  for(Int_t k=0; k<kNDcaVar; k++){
 266    TString histTitle((const char*)fgkResDcaVarTitle[k]);
 267
 268    for(Int_t j=0; j<kNParticles; j++){
 269      for(Int_t i=0; i<kNPtBins; i++){
 270
 271         TString histName((const char*)fgkParticles[j]);
 272         histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 273         TString histEPName((const char*)fgkParticles[j]);
 274         histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 275
 276         if(k==0){
 277           fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 278           fHistDcaXYRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
 279           if(j<(kNParticles-2)){
 280             fHistEPDcaXYRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 281             fHistEPDcaXYRes[j][i]->SetLineColor((const int)fgkColorPart[j]);}
 282         }
 283         if(k==1){
 284           fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 285           fHistDcaZRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
 286           if(j<(kNParticles-2)){
 287             fHistEPDcaZRes[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 288             fHistEPDcaZRes[j][i]->SetLineColor((const int)fgkColorPart[j]); }
 289         }
 290      } // 50 pt bins
 291    } //12 nparticles
 292  } // 2 dca var
 293
 294  //  TList *fResidualList = 0;
 295  residualList->SetOwner();
 296  residualList->SetName("residual");
 297  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 298    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 299      residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
 300      residualList->Add(fHistDcaZRes[iPart][iPtBin]);
 301      if(iPart<(kNParticles-2)){
 302         residualList->Add(fHistEPDcaXYRes[iPart][iPtBin]);
 303         residualList->Add(fHistEPDcaZRes[iPart][iPtBin]);
 304      }
 305    } // loop over pt bins
 306  }  // loop over particles (pos, neg)
 307
 308
 309
 310 }
 311
 312
 313 //________________________________________________________________________
 314 void AliHFEdca::CreateHistogramsPull(TList *pullList){
 315  // define histogram
 316  // 2. pull
 317
 318  const Int_t nBins = 1000;
 319  const Float_t maxXYBin = 20.;
 320  const Float_t maxZBin = 20.;
 321
 322
 323  // for pull -----------------------------------------------------------------------
 324  fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
 325  fHistDcaZPull[kNParticles][kNPtBins]=0x0;
 326
 327  fHistEPDcaXYPull[kNParticles-2][kNPtBins]=0x0;
 328  fHistEPDcaZPull[kNParticles-2][kNPtBins]=0x0;
 329
 330
 331  for(Int_t k=0; k<kNDcaVar; k++){
 332    TString histTitle((const char*)fgkPullDcaVarTitle[k]);
 333
 334    for(Int_t j=0; j<kNParticles; j++){
 335      for(Int_t i=0; i<kNPtBins; i++){
 336
 337         TString histName((const char*)fgkParticles[j]);
 338         histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 339         TString histEPName((const char*)fgkParticles[j]);
 340         histEPName += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 341
 342         if(k==0){
 343           fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
 344           fHistDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 345           if(j<(kNParticles-2))    {
 346             fHistEPDcaXYPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
 347             fHistEPDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);}
 348         }
 349         if(k==1){
 350           fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
 351           fHistDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 352           if(j<(kNParticles-2))    {
 353             fHistEPDcaZPull[j][i] = new TH1F((const char*)histEPName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
 354             fHistEPDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);}
 355         }
 356      } // 50 pt bins
 357    } //6 nparticles
 358  } // 2 dca var
 359
 360  //  TList *fPullList = 0;
 361  pullList->SetOwner();
 362  pullList->SetName("pull");
 363  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 364    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 365      pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
 366      pullList->Add(fHistDcaZPull[iPart][iPtBin]);
 367      if(iPart<(kNParticles-2)){
 368         pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
 369         pullList->Add(fHistDcaZPull[iPart][iPtBin]); }
 370    } // loop over pt bins
 371  }  // loop over particles (pos, neg)
 372
 373 }
 374
 375
 376 //________________________________________________________________________
 377 void AliHFEdca::CreateHistogramsDca(TList *dcaList){
 378  //
 379  // define histograms: MC dca
 380  //
 381
 382  // statistics
 383  fStat = 0x0;
 384  fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
 385  fStat->SetMarkerStyle(20);
 386  fStat->SetMarkerColor(3);
 387  fStat->SetMarkerSize(1);
 388
 389  // for dca
 390  fHistDcaXY[kNParticles][kNPtBins]=0x0;
 391  fHistDcaZ[kNParticles][kNPtBins]=0x0;
 392
 393  fHistEPDcaXY[kNParticles-2][kNPtBins]=0x0;
 394  fHistEPDcaZ[kNParticles-2][kNPtBins]=0x0;
 395
 396  const Int_t nBins = 1000;
 397  const Float_t maxXYBin = 1000.;
 398  const Float_t maxZBin = 1000.;
 399
 400
 401  for(Int_t k=0; k<kNDcaVar; k++){
 402    TString histTitle((const char*)fgkDcaVarTitle[k]);
 403
 404    for(Int_t j=0; j<kNParticles; j++){
 405      for(Int_t i=0; i<kNPtBins; i++){
 406
 407         TString histName((const char*)fgkParticles[j]);
 408         histName += Form("_MCpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 409
 410         TString histNameEP((const char*)fgkParticles[j]);
 411         histNameEP += Form("_ESDpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 412
 413         if(k==0){
 414           fHistDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 415           fHistDcaXY[j][i]->SetLineColor((const int)fgkColorPart[j]);
 416
 417           if(j<(kNParticles-2)){
 418             fHistEPDcaXY[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 419             fHistEPDcaXY[j][i]->SetLineColor((const int)fgkColorPart[j]);}
 420         }
 421         if(k==1){
 422           fHistDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 423           fHistDcaZ[j][i]->SetLineColor((const int)fgkColorPart[j]);
 424           if(j<(kNParticles-2)){
 425             fHistEPDcaZ[j][i] = new TH1F((const char*)histNameEP, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 426             fHistEPDcaZ[j][i]->SetLineColor((const int)fgkColorPart[j]);}
 427         }
 428      } // 50 pt bins
 429    } //12 nparticles
 430  } // 2 dca var
 431
 432    //  TList *fDcaList = 0;
 433  dcaList->SetOwner();
 434  dcaList->SetName("mcDcaDistr");
 435  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 436    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 437      dcaList->Add(fHistDcaXY[iPart][iPtBin]);
 438      dcaList->Add(fHistDcaZ[iPart][iPtBin]);
 439      if(iPart<(kNParticles-2)) {
 440         dcaList->Add(fHistEPDcaXY[iPart][iPtBin]);
 441         dcaList->Add(fHistEPDcaZ[iPart][iPtBin]); }
 442    } // loop over pt bins
 443  }  // loop over particles (pos, neg)
 444
 445  dcaList->Add(fStat);
 446
 447 }
 448
 449 //________________________________________________________________________
 450 void AliHFEdca::CreateHistogramsKfDca(TList *kfDcaList){
 451  //
 452  // define histograms: MC dca
 453  //
 454
 455  // statistics
 456  fStat = 0x0;
 457  fStat = new TH1I("fStatistics", "allStatistics;ID;counts", 7, -3.5, 3.5);
 458  fStat->SetMarkerStyle(20);
 459  fStat->SetMarkerColor(3);
 460  fStat->SetMarkerSize(1);
 461
 462  // for kf dca
 463  fHistKFDcaXY[kNParticles][kNPtBins]=0x0;
 464  fHistKFDcaZ[kNParticles][kNPtBins]=0x0;
 465
 466  const Int_t nBins = 1000;
 467  const Float_t maxXYBin = 1000.;
 468  const Float_t maxZBin = 1000.;
 469
 470
 471  for(Int_t k=0; k<kNDcaVar; k++){
 472    TString histTitle((const char*)fgkDcaVarTitle[k]);
 473
 474    for(Int_t j=0; j<kNParticles; j++){
 475      for(Int_t i=0; i<kNPtBins; i++){
 476         TString histNameKF((const char*)fgkParticles[j]);
 477         histNameKF += Form("_MCpid_KF%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 478
 479         if(k==0){
 480           fHistKFDcaXY[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 481           fHistKFDcaXY[j][i]->SetLineColor((const int)fgkColorPart[j]);
 482         }
 483         if(k==1){
 484           fHistKFDcaZ[j][i] = new TH1F((const char*)histNameKF, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 485           fHistKFDcaZ[j][i]->SetLineColor((const int)fgkColorPart[j]);
 486         }
 487      } // 50 pt bins
 488    } //12 nparticles
 489  } // 2 dca var
 490
 491  kfDcaList->SetOwner();
 492  kfDcaList->SetName("mcKfDcaDistr");
 493  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 494    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 495      kfDcaList->Add(fHistKFDcaXY[iPart][iPtBin]);
 496      kfDcaList->Add(fHistKFDcaZ[iPart][iPtBin]);
 497    } // loop over pt bins
 498  }  // loop over particles (pos, neg)
 499
 500  kfDcaList->Add(fStat);
 501
 502 }
 503
 504
 505 //________________________________________________________________________
 506 void AliHFEdca::CreateHistogramsDataDca(TList *dataDcaList){
 507  //
 508  // define histograms: real Data
 509  //
 510
 511  // for dca
 512  fHistDataDcaXY[kNParticles][kNPtBins]=0x0;
 513  fHistDataDcaZ[kNParticles][kNPtBins]=0x0;
 514
 515  fHistDataWoDcaXY[kNParticles][kNPtBins]=0x0;
 516  fHistDataWoDcaZ[kNParticles][kNPtBins]=0x0;
 517
 518  const Int_t nBins = 1000;
 519  const Float_t maxXYBin = 1000.;
 520  const Float_t maxZBin = 1000.;
 521
 522  for(Int_t k=0; k<kNDcaVar; k++){
 523    TString histTitle((const char*)fgkDcaVarTitle[k]);
 524    for(Int_t j=0; j<kNParticles; j++){
 525      for(Int_t i=0; i<kNPtBins; i++){
 526
 527         TString histName((const char*)fgkParticles[j]);
 528         histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 529
 530         TString histNameWo((const char*)fgkParticles[j]);
 531         histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 532
 533         if(k==0){
 534           fHistDataDcaXY[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 535           fHistDataDcaXY[j][i]->SetLineColor((const int)fgkColorPart[j]);
 536
 537           fHistDataWoDcaXY[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
 538           fHistDataWoDcaXY[j][i]->SetLineColor((const int)fgkColorPart[j]);
 539         }
 540         if(k==1){
 541           fHistDataDcaZ[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 542           fHistDataDcaZ[j][i]->SetLineColor((const int)fgkColorPart[j]);
 543
 544           fHistDataWoDcaZ[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, -maxZBin, maxZBin);
 545           fHistDataWoDcaZ[j][i]->SetLineColor((const int)fgkColorPart[j]);
 546         }
 547      } // 50 pt bins
 548    } //12 nparticles
 549  } // 2 dca var
 550
 551    //  TList *fDcaList = 0;
 552  dataDcaList->SetOwner();
 553  dataDcaList->SetName("dataDcaDistr");
 554  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 555    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 556      dataDcaList->Add(fHistDataDcaXY[iPart][iPtBin]);
 557      dataDcaList->Add(fHistDataDcaZ[iPart][iPtBin]);
 558
 559      dataDcaList->Add(fHistDataWoDcaXY[iPart][iPtBin]);
 560      dataDcaList->Add(fHistDataWoDcaZ[iPart][iPtBin]);
 561    } // loop over pt bins
 562  }  // loop over particles (pos, neg)
 563
 564
 565 }
 566
 567
 568 //________________________________________________________________________
 569 void AliHFEdca::CreateHistogramsDataPull(TList *dataPullList){
 570  // define histogram
 571  // 2. pull
 572
 573  const Int_t nBins = 1000;
 574  const Float_t maxXYBin = 20.;
 575  const Float_t maxZBin = 20.;
 576
 577  // for pull -----------------------------------------------------------------------
 578  fHistDataDcaXYPull[kNParticles][kNPtBins]=0x0;
 579  fHistDataDcaZPull[kNParticles][kNPtBins]=0x0;
 580
 581  fHistDataWoDcaXYPull[kNParticles][kNPtBins]=0x0;
 582  fHistDataWoDcaZPull[kNParticles][kNPtBins]=0x0;
 583
 584
 585  for(Int_t k=0; k<kNDcaVar; k++){
 586    TString histTitle((const char*)fgkPullDataDcaVarTitle[k]);
 587
 588    for(Int_t j=0; j<kNParticles; j++){
 589      for(Int_t i=0; i<kNPtBins; i++){
 590
 591         TString histName((const char*)fgkParticles[j]);
 592         histName += Form("_%s_Data_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 593
 594         TString histNameWo((const char*)fgkParticles[j]);
 595         histNameWo += Form("_%s_Data_wo_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
 596
 597         if(k==0){
 598           fHistDataDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
 599           fHistDataDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 600
 601           fHistDataWoDcaXYPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
 602           fHistDataWoDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 603         }
 604         if(k==1){
 605           fHistDataDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
 606           fHistDataDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 607
 608           fHistDataWoDcaZPull[j][i] = new TH1F((const char*)histNameWo, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
 609           fHistDataWoDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
 610         }
 611      } // 50 pt bins
 612    } //6 nparticles
 613  } // 2 dca var
 614
 615  //  TList *fDataPullList = 0;
 616  dataPullList->SetOwner();
 617  dataPullList->SetName("dataPull");
 618  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 619    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 620      dataPullList->Add(fHistDataDcaXYPull[iPart][iPtBin]);
 621      dataPullList->Add(fHistDataDcaZPull[iPart][iPtBin]);
 622
 623      dataPullList->Add(fHistDataWoDcaXYPull[iPart][iPtBin]);
 624      dataPullList->Add(fHistDataWoDcaZPull[iPart][iPtBin]);
 625    } // loop over pt bins
 626  }  // loop over particles (pos, neg)
 627
 628 }
 629
 630 //________________________________________________________________________
 631 void AliHFEdca::CreateHistogramsVertex(TList *vertexList){
 632  //
 633  // define histograms: vertex
 634  //
 635  // for  vertex
 636
 637  fHistMCvertex[kNVertexVar]=0x0;
 638  fHistESDvertex[kNVertexVar]=0x0;
 639
 640  const Int_t nBins = 1000;
 641  const Float_t minXBin = -0.2e4;
 642  const Float_t maxXBin = 0.2e4;
 643  const Float_t minYBin = -0.5e4;
 644  const Float_t maxYBin = 0.5e4;
 645  const Float_t minZBin = -1.5e5;
 646  const Float_t maxZBin = 1.5e5;
 647
 648  const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
 649  const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
 650
 651  for(Int_t k=0; k<kNVertexVar; k++){
 652    TString histTitle((const char*)fgkVertexVarTitle[k]);
 653    TString histNameMC((const char*)fgkVertexVar[k]);
 654    histNameMC += Form("_MC");
 655    TString histNameESD((const char*)fgkVertexVar[k]);
 656    histNameESD += Form("_ESD");
 657
 658    fHistMCvertex[k] = new TH1F((const char*)histNameMC, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
 659    fHistMCvertex[k]->SetLineColor(k+2);
 660
 661    fHistESDvertex[k] = new TH1F((const char*)histNameESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
 662    fHistESDvertex[k]->SetLineColor(k+2);
 663  } // 3 vertex var
 664
 665  vertexList->SetOwner();
 666  vertexList->SetName("vertexDistr");
 667
 668  for(Int_t k=0; k<kNVertexVar; k++){
 669    vertexList->Add(fHistMCvertex[k]);
 670    vertexList->Add(fHistESDvertex[k]);
 671  }
 672
 673 }
 674
 675
 676
 677 //________________________________________________________________________
 678 void AliHFEdca::CreateHistogramsDataVertex(TList *dataVertexList){
 679  //
 680  // define histograms: vertex
 681  //
 682  // for data vertex
 683
 684  fHistDatavertex[kNVertexVar]=0x0;
 685
 686  const Int_t nBins = 1000;
 687  const Float_t minXBin = -0.2e4;
 688  const Float_t maxXBin = 0.2e4;
 689  const Float_t minYBin = -0.5e4;
 690  const Float_t maxYBin = 0.5e4;
 691  const Float_t minZBin = -1.5e5;
 692  const Float_t maxZBin = 1.5e5;
 693
 694  const Float_t minBin[kNVertexVar] = {minXBin, minYBin, minZBin};
 695  const Float_t maxBin[kNVertexVar] = {maxXBin, maxYBin, maxZBin};
 696
 697  for(Int_t k=0; k<kNVertexVar; k++){
 698    TString histTitle((const char*)fgkVertexVarTitle[k]);
 699    TString histNameDataESD((const char*)fgkVertexVar[k]);
 700    histNameDataESD += Form("_data");
 701
 702    fHistDatavertex[k] = new TH1F((const char*)histNameDataESD, (const char*)histTitle, nBins, minBin[k], maxBin[k]);
 703    fHistDatavertex[k]->SetLineColor(k+2);
 704  } // 3 vertex var
 705
 706    //  TList *fVDaraVertexList = 0;
 707  dataVertexList->SetOwner();
 708  dataVertexList->SetName("dataVertexDistr");
 709
 710  for(Int_t k=0; k<kNVertexVar; k++){
 711    dataVertexList->Add(fHistDatavertex[k]);
 712  }
 713
 714 }
 715
 716 //_______________________________________________________________________________________________
 717 void AliHFEdca::CreateHistogramsPid(TList *mcPidList){
 718  //
 719  // define histograms which fills combined PID
 720  //
 721
 722  const Char_t *mcOResd[2]={"mcPt", "esdPt"};
 723
 724  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 725    TString histTitleMc((const char*)fgkParticles[iPart]);
 726    TString histTitleEsd((const char*)fgkParticles[iPart]);
 727    histTitleMc += Form("_McPid_%s;p_{T} [GeV/c];counts", mcOResd[0]);
 728    histTitleEsd += Form("_EsdPid_%s;p_{T} [GeV/c];counts", mcOResd[1]);
 729
 730    TString histNameMc((const char*)fgkParticles[iPart]);
 731    TString histNameEsd((const char*)fgkParticles[iPart]);
 732    histNameMc+=Form("_McPid_%s", mcOResd[0]);
 733    histNameEsd+=Form("_EsdPid_%s", mcOResd[1]);
 734
 735    fHistMcPid[iPart] = new TH1F(histNameMc, histTitleMc, kNPtBins, fgkPtIntv);
 736    fHistEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
 737  }
 738
 739
 740  mcPidList->SetOwner();
 741  mcPidList->SetName("combinedPid");
 742
 743  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 744    mcPidList->Add(fHistMcPid[iPart]);
 745    mcPidList->Add(fHistEsdPid[iPart]);
 746  }
 747 }
 748
 749
 750 //_______________________________________________________________________________________________
 751 void AliHFEdca::CreateHistogramsDataPid(TList *pidList){
 752  //
 753  // define histograms which fills combined PID: data
 754  //
 755
 756
 757
 758  for(Int_t iPart=0; iPart<kNParticles; iPart++){
 759      TString histTitleEsd((const char*)fgkParticles[iPart]);
 760      histTitleEsd+=Form("_DataEsdPid_esdPt;p_{T} [GeV/c];counts");
 761      TString histNameEsd((const char*)fgkParticles[iPart]);
 762      histNameEsd+=Form("_DataEsdPid");
 763
 764      fHistDataEsdPid[iPart] = new TH1F(histNameEsd, histTitleEsd, kNPtBins, fgkPtIntv);
 765  }
 766
 767
 768  pidList->SetOwner();
 769  pidList->SetName("dataCombinedPid");
 770
 771  for(Int_t iPart=0; iPart<kNParticles; iPart++)
 772    pidList->Add(fHistDataEsdPid[iPart]);
 773
 774
 775 }
 776
 777
 778
 779 //_______________________________________________________________________________________________
 780 void AliHFEdca::FillHistogramsDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
 781 {
 782  // the kDca plugin
 783  // MC vertex
 784  AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
 785  Double_t mcPrimV[3];
 786  mcPrimV[0] = mcPrimVtx->GetX();
 787  mcPrimV[1] = mcPrimVtx->GetY();
 788  mcPrimV[2] = mcPrimVtx->GetZ();
 789
 790  Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
 791
 792 // filling historgams track by track
 793 // obtaining reconstructed dca ------------------------------------------------------------------
 794
 795  Float_t esdpx = track->Px();
 796  Float_t esdpy = track->Py();
 797  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
 798
 799 // obtaining errors of dca ------------------------------------------------------------------
 800  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
 801  Double_t primV[3];
 802  primV[0] = primVtx->GetXv();
 803  primV[1] = primVtx->GetYv();
 804  primV[2] = primVtx->GetZv();
 805
 806  Float_t magneticField = 0;  // initialized as 5kG
 807  magneticField = esdEvent->GetMagneticField();  // in kG
 808
 809  Double_t beampiperadius=3.;
 810  Double_t dz[2];   // error of dca in cm
 811  Double_t covardz[3];
 812
 813  if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return;  // protection
 814  track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz);
 815
 816  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
 817  TParticle *part = mctrack->Particle();
 818
 819  Float_t vx = part->Vx();  // in cm
 820  Float_t vy = part->Vy();  // in cm
 821  Float_t vz = part->Vz();   // in cm
 822
 823  Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
 824
 825  Float_t mcpx = part->Px();
 826  Float_t mcpy = part->Py();
 827  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
 828
 829  Int_t pdg = part->GetPdgCode();
 830  Int_t esdPid = GetCombinedPid(track);
 831
 832  Int_t charge = 1;
 833  if(pdg==kPDGelectron || pdg==kPDGmuon
 834     || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
 835
 836  // calculate mcDca ------------------------------------------------------------------
 837  const Float_t conv[2] = {1.783/1.6, 2.99792458};
 838  Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
 839
 840  Float_t nx = esdpx/mcpt;
 841  Float_t ny = esdpy/mcpt;
 842
 843  Float_t radius;
 844  radius = TMath::Abs(radiusMc);
 845
 846  Double_t dxy = vxy - mcVtxXY;   // in cm
 847  Double_t dvx = vx - mcPrimV[0]; // in cm
 848  Double_t dvy = vy - mcPrimV[1]; // in cm
 849
 850  Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
 851
 852  Double_t mcDca[2] = {mcDcaXY*100, vz};  // in cm
 853  Double_t residual[2] = {0, 0};
 854  Double_t pull[2] = {0, 0};
 855  Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
 856  for(Int_t i=0; i<2; i++){
 857    residual[i] = dz[i] - mcDca[i]; // in centimeters
 858    if(error[i]!=0)pull[i] = residual[i]/error[i];   // unitless
 859  }
 860
 861
 862  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
 863    // identified ones
 864    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 865      if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
 866         if(pdg==fgkPdgParticle[iPart]) {
 867           fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
 868           fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
 869           fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
 870           fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
 871           fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
 872           fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
 873         }  // mc pdg
 874
 875         if(esdPid==fgkPdgParticle[iPart]) {
 876           fHistEPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
 877           fHistEPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
 878           fHistEPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
 879           fHistEPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
 880           fHistEPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
 881           fHistEPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
 882         }  // esd pid
 883
 884      } // pt range
 885
 886      else
 887         continue;
 888    }  // pt loop
 889  } // particle id loop
 890
 891  // for charged particles: no pid
 892  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
 893    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
 894      Int_t iPart = 10;
 895      if(charge>0) iPart = 11;
 896      fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
 897      fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
 898      fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
 899      fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
 900      fHistDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
 901      fHistDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
 902    }
 903    else
 904      continue;
 905  } // pt
 906
 907 }
 908
 909 //_______________________________________________________________________________________________
 910 void AliHFEdca::FillHistogramsKfDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
 911 {
 912  // the kKfDca plugin
 913
 914 // filling historgams track by track
 915
 916 // obtaining reconstructed dca ------------------------------------------------------------------
 917  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
 918    Float_t magneticField = 0;  // initialized as 5kG
 919  magneticField = esdEvent->GetMagneticField();  // in kG
 920
 921  Float_t esdpx = track->Px();
 922  Float_t esdpy = track->Py();
 923  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
 924
 925  Int_t charge = track->Charge();
 926
 927  Double_t beampiperadius=3.;
 928  Double_t dz[2];   // error of dca in cm
 929  Double_t covardz[3];
 930  if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
 931  track->PropagateToDCA(primVtx,magneticField,  beampiperadius, dz, covardz);
 932
 933  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
 934  TParticle *part = mctrack->Particle();
 935  Int_t pdg = part->GetPdgCode();
 936
 937  // calculate dca using AliKFParticle class------------------------------------------------------------------
 938  Double_t kfdz[3] = {0, 0, 0};
 939  Double_t kfdzwith[3] = {0, 0, 0};
 940
 941  Int_t trkID = track->GetID();
 942
 943  AliKFParticle::SetField(magneticField);
 944  AliKFParticle kfParticle(*track, pdg);
 945
 946  // prepare kfprimary vertex
 947  AliKFVertex kfESDprimary;
 948  // Reconstruct Primary Vertex (with ESD tracks)
 949  Int_t n=primVtx->GetNIndices();
 950  if (n>0 && primVtx->GetStatus()){
 951    kfESDprimary = AliKFVertex(*primVtx);
 952
 953    Double_t dcaXYWithTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
 954    Double_t dcaWithTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
 955    Double_t dcaZWithTrk = 0;
 956    if(TMath::Abs(dcaWithTrk)>=TMath::Abs(dcaXYWithTrk))
 957      dcaZWithTrk =TMath::Sqrt(dcaWithTrk*dcaWithTrk-dcaXYWithTrk*dcaXYWithTrk)*((dz[1]*-1<=0)?1:-1);
 958    kfdzwith[0] = dcaXYWithTrk;
 959    kfdzwith[1] = dcaZWithTrk;
 960    kfdzwith[2] = dcaWithTrk;  // with current track
 961
 962    Double_t dcaXYWoTrk = 0;
 963    Double_t dcaZWoTrk = 0;
 964    Double_t dcaWoTrk = 0;
 965
 966    UShort_t *priIndex = primVtx->GetIndices();
 967
 968    for (Int_t i=0;i<n;i++){
 969
 970      Int_t idx = Int_t(priIndex[i]);
 971      if (idx == trkID){
 972         kfESDprimary -= kfParticle;
 973         dcaXYWoTrk = kfParticle.GetDistanceFromVertexXY(kfESDprimary);
 974         dcaWoTrk = kfParticle.GetDistanceFromVertex(kfESDprimary);
 975         if((dcaWoTrk-dcaXYWoTrk)>=0)
 976           dcaZWoTrk = TMath::Abs(dcaWoTrk*dcaWoTrk - dcaXYWoTrk*dcaXYWoTrk)*((dz[1]*-1<=0)?1:-1);
 977      }  // remove current track from this calculation
 978    }  // loop over all primary vertex contributors
 979
 980
 981    kfdz[0] = dcaXYWoTrk;
 982    kfdz[1] = dcaZWoTrk;
 983    kfdz[2] = dcaWoTrk;
 984
 985  }  // only if n contributor > 0 and primVtx constructed
 986
 987  fStat->Fill(0);
 988
 989  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
 990  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
 991  if(kfdzwith[0]==0 && dz[0]!=0) fStat->Fill(3);  // 0 from KF particle (with current track)
 992
 993  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
 994  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
 995  if(kfdz[0]==0 && dz[0]!=0) fStat->Fill(-3);  // 0 from KF particle (without current track)
 996
 997  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
 998    // identified ones
 999    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1000      if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1001         if(pdg==fgkPdgParticle[iPart]) {
1002           fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1.0e4);
1003           fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1.0e4);
1004         }  // mc pdg
1005      } // pt range
1006
1007      else
1008         continue;
1009    }  // pt loop
1010  } // particle id loop
1011
1012  // for charged particles: no pid
1013  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1014    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1015      Int_t iPart = 10;
1016      if(charge>0) iPart = 11;
1017      fHistKFDcaXY[iPart][iPtBin]->Fill(kfdzwith[0]*1e4);
1018      fHistKFDcaZ[iPart][iPtBin]->Fill(kfdzwith[1]*1e4);
1019
1020    }
1021    else
1022      continue;
1023  } // pt
1024
1025 }  // KF dca
1026
1027
1028 //_______________________________________________________________________________________________
1029 void AliHFEdca::FillHistogramsVtx(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
1030 {
1031
1032  // MC vertex
1033  AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1034  Double_t mcPrimV[3];
1035  mcPrimV[0] = mcPrimVtx->GetX();
1036  mcPrimV[1] = mcPrimVtx->GetY();
1037  mcPrimV[2] = mcPrimVtx->GetZ();
1038
1039 // obtaining errors of dca ------------------------------------------------------------------
1040  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1041  Double_t primV[3];
1042  primV[0] = primVtx->GetXv();
1043  primV[1] = primVtx->GetYv();
1044  primV[2] = primVtx->GetZv();
1045
1046  for(Int_t i=0; i<kNVertexVar; i++){
1047    fHistMCvertex[i]->Fill(mcPrimV[i]*1.0e4);
1048    fHistESDvertex[i]->Fill(primV[i]*1.0e4);
1049  }
1050
1051 }
1052
1053 //_______________________________________________________________________________________________
1054 void AliHFEdca::FillHistogramsPid(AliESDtrack * const track, AliMCEvent * const mcEvent)
1055 {
1056
1057
1058 // filling historgams track by track
1059  Float_t esdpx = track->Px();
1060  Float_t esdpy = track->Py();
1061  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1062
1063  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1064  TParticle *part = mctrack->Particle();
1065
1066  Float_t mcpx = part->Px();
1067  Float_t mcpy = part->Py();
1068  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1069
1070  Int_t pdg = part->GetPdgCode();
1071  Int_t esdPid = GetCombinedPid(track);
1072
1073
1074  Double_t ptMom[2] = {mcpt, esdpt};
1075  // for combined PID
1076  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1077    if(pdg==fgkPdgParticle[iPart])                // pid all by MC
1078      fHistMcPid[iPart]->Fill(ptMom[0]);
1079
1080    if(esdPid==fgkPdgParticle[iPart])             // pid all by combined pid
1081      fHistEsdPid[iPart]->Fill(ptMom[1]);
1082  } // loop over particles
1083
1084  // for charged
1085  if(pdg==kPDGelectron || pdg==kPDGmuon || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton)
1086    fHistMcPid[10]->Fill(ptMom[0]);
1087  if(pdg==-kPDGelectron || pdg==-kPDGmuon || pdg==kPDGpion || pdg==kPDGkaon || pdg==kPDGproton)
1088    fHistMcPid[11]->Fill(ptMom[0]);
1089  if(esdPid==kPDGelectron || esdPid==kPDGmuon || esdPid==-kPDGpion || esdPid==-kPDGkaon || esdPid==-kPDGproton)
1090    fHistEsdPid[10]->Fill(ptMom[1]);
1091  if(esdPid==-kPDGelectron || esdPid==-kPDGmuon || esdPid==kPDGpion || esdPid==kPDGkaon || esdPid==kPDGproton)
1092    fHistEsdPid[11]->Fill(ptMom[1]);
1093 }
1094
1095
1096 ////_______________________________________________________________________________________________
1097 void AliHFEdca::FillHistogramsDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliESDVertex * const vtxESDSkip)
1098 {
1099 // filling historgams track by track
1100 // obtaining reconstructed dca --------------------------------------------------------------
1101
1102  Float_t esdpx = track->Px();
1103  Float_t esdpy = track->Py();
1104  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1105  Int_t charge = track->Charge();
1106
1107 // obtaining errors of dca ------------------------------------------------------------------
1108  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1109  Double_t primV[3];
1110  primV[0] = primVtx->GetXv();
1111  primV[1] = primVtx->GetYv();
1112  primV[2] = primVtx->GetZv();
1113
1114
1115  Float_t magneticField = 0;  // initialized as 5kG
1116  magneticField = esdEvent->GetMagneticField();  // in kG
1117
1118  Double_t beampiperadius=3.;
1119  Double_t dz[2];   // error of dca in cm
1120  Double_t covardz[3];
1121
1122  if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return;  // protection
1123  track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz);
1124
1125  Double_t pull[2] = {0, 0};
1126  Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1127  for(Int_t i=0; i<2; i++){
1128    if(error[i]!=0)pull[i] = dz[i]/error[i];   // unitless
1129  }
1130
1131  // get dca when current track is not included
1132
1133  Double_t dzwo[2], covardzwo[3];
1134  Double_t pullwo[2] = {0, 0};
1135  if(!track->PropagateToDCA(vtxESDSkip, magneticField, beampiperadius, dzwo, covardzwo)) return;   // protection
1136  track->PropagateToDCA(vtxESDSkip, magneticField, beampiperadius, dzwo, covardzwo);
1137
1138  Double_t errorwo[2] ={TMath::Sqrt(TMath::Abs(covardzwo[0])), TMath::Sqrt(TMath::Abs(covardzwo[2]))};
1139  for(Int_t i=0; i<2; i++){
1140    if(errorwo[i]!=0) pullwo[i] = dzwo[i]/errorwo[i];   // unitless
1141  }
1142
1143  // do pid
1144  Int_t esdPid = GetCombinedPid(track);
1145
1146  for(Int_t iPart=0; iPart<(kNParticles-2); iPart++){
1147    // identified ones
1148    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1149      if(esdPid==fgkPdgParticle[iPart] && (esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1])) {
1150         fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1151         fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1152         fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1153         fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1154         // w/o current track
1155         fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1156         fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1157         fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1158         fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1159      }
1160      else
1161         continue;
1162    }
1163  }
1164
1165  // for charged particles
1166  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1167    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1168      Int_t iPart = 10;
1169      if(charge>0) iPart = 11;
1170      fHistDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1171      fHistDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1172      fHistDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1173      fHistDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1174      // without current track
1175      fHistDataWoDcaXY[iPart][iPtBin]->Fill(dzwo[0]*1e4);
1176      fHistDataWoDcaZ[iPart][iPtBin]->Fill(dzwo[1]*1e4);
1177      fHistDataWoDcaXYPull[iPart][iPtBin]->Fill(pullwo[0]);
1178      fHistDataWoDcaZPull[iPart][iPtBin]->Fill(pullwo[1]);
1179
1180    }
1181    else
1182      continue;
1183  }
1184
1185 }
1186
1187 //_______________________________________________________________________________________________
1188 void AliHFEdca::FillHistogramsDataVtx(AliESDEvent * const esdEvent)
1189 {
1190
1191
1192 // obtaining errors of dca ------------------------------------------------------------------
1193  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1194  Double_t primV[3];
1195  primV[0] = primVtx->GetXv();
1196  primV[1] = primVtx->GetYv();
1197  primV[2] = primVtx->GetZv();
1198
1199  // require events with at least 3 contributors for primary vertex construction
1200  Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1201  if(nEsdPrimVtxCtrb<1) return; // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1202  for(Int_t i=0; i<kNVertexVar; i++)
1203    fHistDatavertex[i]->Fill(primV[i]*1e4);
1204
1205 }
1206
1207
1208 ////_______________________________________________________________________________________________
1209 void AliHFEdca::FillHistogramsDataPid(AliESDtrack * const track)
1210 {
1211 // filling historgams track by track
1212 // obtaining reconstructed dca --------------------------------------------------------------
1213
1214  Float_t esdpx = track->Px();
1215  Float_t esdpy = track->Py();
1216  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1217  Int_t charge = track->Charge();
1218
1219  Int_t esdPid = GetCombinedPid(track);
1220
1221  for(Int_t iPart=0; iPart<kNParticles; iPart++){
1222    if(iPart<(kNParticles-2)){
1223      if(esdPid==fgkPdgParticle[iPart])  fHistDataEsdPid[iPart]->Fill(esdpt);
1224    }     // for identified
1225
1226    else {
1227      if(charge<0) fHistDataEsdPid[10]->Fill(esdpt);
1228      if(charge>0) fHistDataEsdPid[11]->Fill(esdpt);
1229    }
1230  }
1231 }
1232
1233 //_________________________________________________________________________________________________
1234 void AliHFEdca::ApplyExtraCuts(AliESDEvent * const esdEvent, Int_t nMinPrimVtxContributor)
1235 {
1236
1237  //
1238  // only one extra cut, number of contributors to each primary vertex
1239  //
1240
1241  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1242  Int_t nEsdPrimVtxCtrb = primVtx->GetNContributors();
1243  if(nEsdPrimVtxCtrb<nMinPrimVtxContributor) return;
1244  // for pass 1, no diomond constrain, each event has at least 1 contributor to Vtx
1245
1246 }
1247
1248 //_____________________________________________________
1249 Int_t AliHFEdca::GetCombinedPid(AliESDtrack *const track)
1250 {
1251
1252  //combined detector pid
1253  Double_t prob[AliPID::kSPECIES];
1254  track->GetESDpid(prob);
1255  Double_t priors[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
1256
1257  Int_t charge = track->Charge();
1258  Int_t esdPid = -1;
1259
1260  AliPID pid;
1261  pid.SetPriors(priors);
1262  pid.SetProbabilities(prob);
1263
1264  // identify particle as the most probable
1265
1266  Double_t pelectron = pid.GetProbability(AliPID::kElectron);
1267  if(pelectron > pid.GetProbability(AliPID::kMuon) &&
1268     pelectron > pid.GetProbability(AliPID::kPion) &&
1269     pelectron > pid.GetProbability(AliPID::kKaon) &&
1270     pelectron > pid.GetProbability(AliPID::kProton) )  esdPid = -kPDGelectron;
1271
1272  Double_t pmuon = pid.GetProbability(AliPID::kMuon);
1273  if(pmuon > pid.GetProbability(AliPID::kElectron) &&
1274     pmuon > pid.GetProbability(AliPID::kPion) &&
1275     pmuon > pid.GetProbability(AliPID::kKaon) &&
1276     pmuon > pid.GetProbability(AliPID::kProton) )  esdPid = -kPDGmuon;
1277
1278  Double_t ppion = pid.GetProbability(AliPID::kPion);
1279  if(ppion > pid.GetProbability(AliPID::kElectron) &&
1280     ppion > pid.GetProbability(AliPID::kMuon) &&
1281     ppion > pid.GetProbability(AliPID::kKaon) &&
1282     ppion > pid.GetProbability(AliPID::kProton) )  esdPid = kPDGpion;
1283
1284  Double_t pkaon = pid.GetProbability(AliPID::kKaon);
1285  if(pkaon > pid.GetProbability(AliPID::kElectron) &&
1286     pkaon > pid.GetProbability(AliPID::kMuon) &&
1287     pkaon > pid.GetProbability(AliPID::kPion) &&
1288     pkaon > pid.GetProbability(AliPID::kProton) )  esdPid = kPDGkaon;
1289
1290  Double_t pproton = pid.GetProbability(AliPID::kProton);
1291  if(pproton > pid.GetProbability(AliPID::kElectron) &&
1292     pproton > pid.GetProbability(AliPID::kMuon) &&
1293     pproton > pid.GetProbability(AliPID::kPion) &&
1294     pproton > pid.GetProbability(AliPID::kKaon) )  esdPid = kPDGproton;
1295
1296
1297  return charge*esdPid;
1298
1299 }
1300
1301
1302 // for the HFE pid
1303
1304 //________________________________________________________________________
1305 void AliHFEdca::CreateHistogramsHfeDca(TList *hfeDcaList){
1306  //
1307  // define histograms: hfe pid electrons in MC
1308  //
1309
1310  const Int_t nBinsDca = 1000;
1311  const Float_t maxXYBinDca = 1000.;
1312  const Float_t maxZBinDca = 1000.;
1313
1314  const Int_t nBinsPull = 1000;
1315  const Float_t maxXYBinPull = 20.;
1316  const Float_t maxZBinPull = 20.;
1317
1318  const Char_t *mcOResd[2]={"mcPt", "esdPt"};
1319
1320  fHistHPDcaXY[2][kNPtBins]=0x0;
1321  fHistHPDcaZ[2][kNPtBins]=0x0;
1322  fHistHPDcaXYRes[2][kNPtBins]=0x0;
1323  fHistHPDcaZRes[2][kNPtBins]=0x0;
1324  fHistHPDcaXYPull[2][kNPtBins]=0x0;
1325  fHistHPDcaZPull[2][kNPtBins]=0x0;
1326
1327
1328  for(Int_t k=0; k<kNDcaVar; k++){
1329
1330    TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1331    TString histTitleRes((const char*)fgkPullDcaVarTitle[k]);
1332    TString histTitlePull((const char*)fgkResDcaVarTitle[k]);
1333
1334    for(Int_t iPart=0; iPart<2; iPart++){
1335      for(Int_t i=0; i<kNPtBins; i++){
1336         TString histHPName((const char*)fgkParticles[iPart*5]);
1337         histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1338         TString histHPNameRes((const char*)fgkParticles[iPart*5]);
1339         histHPNameRes += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkResDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1340         TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1341         histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1342
1343         if(k==0){
1344           fHistHPDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1345           fHistHPDcaXY[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1346           fHistHPDcaXYRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1347           fHistHPDcaXYRes[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1348           fHistHPDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1349           fHistHPDcaXYPull[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1350         }
1351
1352         if(k==1){
1353           fHistHPDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1354           fHistHPDcaZ[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1355           fHistHPDcaZRes[iPart][i] = new TH1F((const char*)histHPNameRes, (const char*)histTitleRes, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1356           fHistHPDcaZRes[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1357           fHistHPDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxZBinPull, 1+maxZBinPull);
1358           fHistHPDcaZPull[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1359         }
1360      } // 50 pt bins
1361    } //2 nparticles
1362  } // 2 dca var
1363
1364  fHistHfePid[2][2] = 0x0; //!  HFE pid
1365  for(Int_t id=0; id<2; id++){
1366    for(Int_t iPart=0; iPart<2; iPart++){
1367      TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1368      histTitleHfe+=Form("_MC_HfePid_esdPt;p_{T} [GeV/c];counts");
1369      TString histNameHfe((const char*)fgkParticles[iPart*5]);
1370      histNameHfe+=Form("_MC_HfePid_%s", mcOResd[id]);
1371      fHistHfePid[id][iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1372    }
1373  }
1374
1375  hfeDcaList->SetOwner();
1376  hfeDcaList->SetName("hfeDca");
1377  for(Int_t iPart=0; iPart<2; iPart++){
1378    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1379      hfeDcaList->Add(fHistHPDcaXY[iPart][iPtBin]);
1380      hfeDcaList->Add(fHistHPDcaZ[iPart][iPtBin]);
1381      hfeDcaList->Add(fHistHPDcaXYRes[iPart][iPtBin]);
1382      hfeDcaList->Add(fHistHPDcaZRes[iPart][iPtBin]);
1383      hfeDcaList->Add(fHistHPDcaXYPull[iPart][iPtBin]);
1384      hfeDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1385    }
1386    for(Int_t id=0; id<2; id++)
1387      hfeDcaList->Add(fHistHfePid[id][iPart]);
1388  }
1389
1390 }
1391
1392
1393 //________________________________________________________________________
1394 void AliHFEdca::CreateHistogramsHfeDataDca(TList *hfeDataDcaList){
1395  //
1396  // define histograms: hfe pid electrons in data
1397  //
1398
1399  const Int_t nBinsDca = 1000;
1400  const Float_t maxXYBinDca = 1000.;
1401  const Float_t maxZBinDca = 1000.;
1402
1403  const Int_t nBinsPull = 1000;
1404  const Float_t maxXYBinPull = 20.;
1405  const Float_t maxZBinPull = 20.;
1406
1407
1408  fHistHPDataDcaXY[2][kNPtBins]=0x0;
1409  fHistHPDataDcaZ[2][kNPtBins]=0x0;
1410  fHistHPDataDcaXYPull[2][kNPtBins]=0x0;
1411  fHistHPDataDcaZPull[2][kNPtBins]=0x0;
1412
1413  for(Int_t k=0; k<kNDcaVar; k++){
1414    TString histTitleDca((const char*)fgkDcaVarTitle[k]);
1415    TString histTitlePull((const char*)fgkPullDcaVarTitle[k]);
1416    for(Int_t iPart=0; iPart<2; iPart++){
1417      for(Int_t i=0; i<kNPtBins; i++){
1418         TString histHPName((const char*)fgkParticles[iPart*5]);
1419         histHPName += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1420         TString histHPNamePull((const char*)fgkParticles[iPart*5]);
1421         histHPNamePull += Form("_HFEpid_%s_pT-%.2f-%.2f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
1422
1423         if(k==0){
1424           fHistHPDataDcaXY[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxXYBinDca, 1+maxXYBinDca);
1425           fHistHPDataDcaXY[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1426           fHistHPDataDcaXYPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsPull, 1-maxXYBinPull, 1+maxXYBinPull);
1427           fHistHPDataDcaXYPull[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1428         }
1429
1430         if(k==1){
1431           fHistHPDataDcaZ[iPart][i] = new TH1F((const char*)histHPName, (const char*)histTitleDca, nBinsDca, 1-maxZBinDca, 1+maxZBinDca);
1432           fHistHPDataDcaZ[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1433           fHistHPDataDcaZPull[iPart][i] = new TH1F((const char*)histHPNamePull, (const char*)histTitlePull, nBinsDca, 1-maxZBinPull, 1+maxZBinPull);
1434           fHistHPDataDcaZPull[iPart][i]->SetLineColor((const int)fgkColorPart[iPart*5]);
1435         }
1436
1437      } // 50 pt bins
1438    } // 2 particle type
1439  } // 2 dca var
1440
1441  fHistDataHfePid[2] = 0x0; //!  HFE pid
1442  for(Int_t iPart=0; iPart<2; iPart++){
1443    TString histTitleHfe((const char*)fgkParticles[iPart*5]);
1444    histTitleHfe+=Form("_Data_HfePid_esdPt;p_{T} [GeV/c];counts");
1445    TString histNameHfe((const char*)fgkParticles[iPart*5]);
1446    histNameHfe+=Form("_Data_HfePid");
1447    fHistDataHfePid[iPart] = new TH1F(histNameHfe, histTitleHfe, kNPtBins, fgkPtIntv);
1448  }
1449
1450
1451  hfeDataDcaList->SetOwner();
1452  hfeDataDcaList->SetName("hfeDataDca");
1453  for(Int_t iPart=0; iPart<2; iPart++){
1454    for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1455      hfeDataDcaList->Add(fHistHPDataDcaXY[iPart][iPtBin]);
1456      hfeDataDcaList->Add(fHistHPDataDcaZ[iPart][iPtBin]);
1457      hfeDataDcaList->Add(fHistHPDataDcaXYPull[iPart][iPtBin]);
1458      hfeDataDcaList->Add(fHistHPDcaZPull[iPart][iPtBin]);
1459
1460      hfeDataDcaList->Add(fHistDataHfePid[iPart]);
1461    }
1462  }
1463
1464
1465 }
1466
1467
1468 //_______________________________________________________________________________________________
1469 void AliHFEdca::FillHistogramsHfeDca(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
1470 {
1471  // the kHFEpid plugin
1472
1473  AliMCVertex *mcPrimVtx = (AliMCVertex *)mcEvent->GetPrimaryVertex();
1474  Double_t mcPrimV[3];
1475  mcPrimV[0] = mcPrimVtx->GetX();
1476  mcPrimV[1] = mcPrimVtx->GetY();
1477  mcPrimV[2] = mcPrimVtx->GetZ();
1478
1479  Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
1480
1481 // filling historgams track by track
1482 // obtaining reconstructed dca ------------------------------------------------------------------
1483  Float_t esdpx = track->Px();
1484  Float_t esdpy = track->Py();
1485  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1486
1487 // obtaining errors of dca ------------------------------------------------------------------
1488  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1489  Double_t primV[3];
1490  primV[0] = primVtx->GetXv();
1491  primV[1] = primVtx->GetYv();
1492  primV[2] = primVtx->GetZv();
1493
1494  Float_t magneticField = 0;  // initialized as 5kG
1495  magneticField = esdEvent->GetMagneticField();  // in kG
1496
1497  Double_t beampiperadius=3.;
1498  Double_t dz[2];   // error of dca in cm
1499  Double_t covardz[3];
1500  if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1501  track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz);
1502
1503  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
1504  TParticle *part = mctrack->Particle();
1505
1506  Float_t vx = part->Vx();  // in cm
1507  Float_t vy = part->Vy();  // in cm
1508  Float_t vz = part->Vz();   // in cm
1509
1510  Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
1511
1512  Float_t mcpx = part->Px();
1513  Float_t mcpy = part->Py();
1514  Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
1515
1516  Int_t pdg = part->GetPdgCode();
1517
1518  Int_t charge = 1;
1519  if(pdg==kPDGelectron || pdg==kPDGmuon
1520     || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
1521
1522  // calculate mcDca ------------------------------------------------------------------
1523  const Float_t conv[2] = {1.783/1.6, 2.99792458};
1524  Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
1525
1526  Float_t nx = esdpx/mcpt;
1527  Float_t ny = esdpy/mcpt;
1528
1529  Float_t radius;
1530  radius = TMath::Abs(radiusMc);
1531
1532  Double_t dxy = vxy - mcVtxXY;   // in cm
1533  Double_t dvx = vx - mcPrimV[0]; // in cm
1534  Double_t dvy = vy - mcPrimV[1]; // in cm
1535
1536  Float_t mcDcaXY = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
1537
1538  Double_t mcDca[2] = {mcDcaXY*100, vz};  // in cm
1539  Double_t residual[2] = {0, 0};
1540  Double_t pull[2] = {0, 0};
1541  Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1542  for(Int_t i=0; i<2; i++){
1543    residual[i] = dz[i] - mcDca[i]; // in centimeters
1544    if(error[i]!=0)pull[i] = residual[i]/error[i];   // unitless
1545  }
1546
1547  Int_t iPart = -1;
1548  if(track->Charge()<0) iPart = 0;  // electron
1549  if(track->Charge()>0) iPart = 1;  // positron
1550  if(track->Charge()==0) {
1551    printf("this is not an electron! Check HFEpid method");
1552    return;
1553  }
1554  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1555    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
1556      fHistHPDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4);
1557      fHistHPDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4);
1558      fHistHPDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1559      fHistHPDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1560      fHistHPDcaXY[iPart][iPtBin]->Fill(dz[0]*1.0e4);
1561      fHistHPDcaZ[iPart][iPtBin]->Fill(dz[1]*1.0e4);
1562
1563    } // pt range
1564
1565    else
1566      continue;
1567  }  // pt loop
1568
1569  fHistHfePid[iPart][0]->Fill(esdpt);
1570  fHistHfePid[iPart][1]->Fill(mcpt);
1571
1572 }
1573
1574
1575 //_______________________________________________________________________________________________
1576 void AliHFEdca::FillHistogramsHfeDataDca(AliESDEvent * const esdEvent, AliESDtrack * const track)
1577 {
1578 // filling historgams track by track
1579 // obtaining reconstructed dca --------------------------------------------------------------
1580
1581  Float_t esdpx = track->Px();
1582  Float_t esdpy = track->Py();
1583  Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
1584  Int_t charge = track->Charge();
1585
1586 // obtaining errors of dca ------------------------------------------------------------------
1587  const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
1588  Double_t primV[3];
1589  primV[0] = primVtx->GetXv();
1590  primV[1] = primVtx->GetYv();
1591  primV[2] = primVtx->GetZv();
1592
1593  Float_t magneticField = 0;  // initialized as 5kG
1594  magneticField = esdEvent->GetMagneticField();  // in kG
1595  Double_t beampiperadius=3.;
1596  Double_t dz[2];   // error of dca in cm
1597  Double_t covardz[3];
1598  if(!track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz)) return; // protection
1599  track->PropagateToDCA(primVtx,magneticField, beampiperadius, dz, covardz);
1600
1601  Double_t pull[2] = {0, 0};
1602  Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
1603  for(Int_t i=0; i<2; i++){
1604    if(error[i]!=0)pull[i] = dz[i]/error[i];   // unitless
1605  }
1606
1607  Int_t iPart = -1;
1608  if(charge<0) iPart = 0;  // electron
1609  if(charge>0) iPart = 1;  // positron
1610  if(charge==0) {
1611    printf("this is not an electron! Check HFEpid method\n");
1612    return;
1613  }
1614
1615  for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
1616    if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]) {
1617      fHistHPDataDcaXY[iPart][iPtBin]->Fill(dz[0]*1e4);
1618      fHistHPDataDcaZ[iPart][iPtBin]->Fill(dz[1]*1e4);
1619      fHistHPDataDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
1620      fHistHPDataDcaZPull[iPart][iPtBin]->Fill(pull[1]);
1621
1622    }
1623    else
1624      continue;
1625  }
1626
1627  fHistDataHfePid[iPart]->Fill(esdpt);
1628
1629 }
1630