PWGHF/hfe/AliHFEdca.cxx

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