PWG3/hfe/AliHFEdisplacedElectrons.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 electrons from beauty study
  17 // Counting electrons from beauty
  18 // by DCA cuts, background subtraction
  19 //
  20 // Authors:
  21 //   Hongyan Yang <hongyan@physi.uni-heidelberg.de>
  22 //
  23
  24 #include "TMath.h"
  25 #include "TList.h"
  26 #include "AliLog.h"
  27
  28 #include <TParticle.h>
  29 #include <TDatabasePDG.h>
  30
  31 #include "THnSparse.h"
  32 #include "AliMCEvent.h"
  33 #include "AliESDEvent.h"
  34 #include "AliMCParticle.h"
  35 #include "AliESDtrack.h"
  36
  37 #include "AliStack.h"
  38
  39 #include "AliHFEdisplacedElectrons.h"
  40
  41 ClassImp(AliHFEdisplacedElectrons)
  42
  43 //__________________________________________________________
  44 const Float_t AliHFEdisplacedElectrons::fgkDcaMinPtIntv[13] = {
  45   //
  46   // define DCA low limits for single electrons cut in each Pt bin
  47   // preliminary numbers
  48   // these numbers should be replaced by the best numbers determined by
  49   // cut efficiency study: signal/background  (not used right now)
  50   //
  51   450, // 0.0 - 0.5
  52   450, // 0.5 - 1.0
  53   450, // 1.0 - 1.5
  54   500, // 1.5 - 2.0
  55   400, // 2.0 - 2.5
  56   300, // 2.5 - 3.0
  57   300, // 3.0 - 4.0
  58   300, // 4.0 - 5.0
  59   200, // 5.0 - 7.0
  60   150, // 7.0 - 9.0
  61   150, // 9.0 - 12.0
  62   100, //12.0 - 16.0
  63   50}; //16.0 - 20.0
  64 //__________________________________________________________
  65 const Float_t AliHFEdisplacedElectrons::fgkPtIntv[14] = {
  66   //
  67   // define pT bins for spectra of single electrons
  68   //
  69   0.0,0.5,1.0,1.5,2.0,2.5,3.0,4.0,5.0,7.0,9.0,12.0,16.0,20.0};
  70
  71 //__________________________________________________________
  72 const Char_t *AliHFEdisplacedElectrons::fgkKineVar[3] = {
  73   "y", "phi", "pt"};
  74
  75 //__________________________________________________________
  76 const Char_t *AliHFEdisplacedElectrons::fgkKineVarTitle[3] ={
  77   "rapidity;y;dN/dy",   "azimuthal;#phi;dN/d#phi",   "transverse momentum;p_{T};dN/dp_{T}",
  78 };
  79
  80 //__________________________________________________________
  81 AliHFEdisplacedElectrons::AliHFEdisplacedElectrons():
  82   fDebugLevel(0)
  83   , fESD(0x0)
  84   , fMC(0x0)
  85   , fTHnSparseDcaMcPionInfo(NULL)
  86   , fTHnSparseDcaMcEleInfo(NULL)
  87   , fTHnSparseDcaDataEleInfo(NULL)
  88   , fOutputList(0x0)
  89 {
  90   //
  91   // default constructor
  92   //
  93
  94 }
  95
  96 //__________________________________________________________
  97 AliHFEdisplacedElectrons::AliHFEdisplacedElectrons(const AliHFEdisplacedElectrons &displacedElectrons):
  98   TObject(displacedElectrons)
  99   , fDebugLevel(0)
 100   , fESD(0x0)
 101   , fMC(0x0)
 102   , fTHnSparseDcaMcPionInfo(NULL)
 103   , fTHnSparseDcaMcEleInfo(NULL)
 104   , fTHnSparseDcaDataEleInfo(NULL)
 105   , fOutputList(0x0)
 106
 107 {
 108   //
 109   // copy constructor
 110   //
 111 }
 112
 113
 114 //__________________________________________________________
 115 AliHFEdisplacedElectrons&AliHFEdisplacedElectrons::operator=(const AliHFEdisplacedElectrons &)
 116 {
 117   //
 118   // Assignment operator
 119   //
 120   Printf("Not yet implemented.");
 121   return *this;
 122 }
 123
 124 //__________________________________________________________
 125 AliHFEdisplacedElectrons::~AliHFEdisplacedElectrons()
 126 {
 127   //
 128   // default constructor
 129   //
 130
 131   if(fTHnSparseDcaMcPionInfo)
 132     delete fTHnSparseDcaMcPionInfo;
 133   if(fTHnSparseDcaMcEleInfo)
 134     delete fTHnSparseDcaMcEleInfo;
 135   if(fTHnSparseDcaDataEleInfo)
 136     delete fTHnSparseDcaDataEleInfo;
 137
 138   if(fOutputList){
 139     fOutputList->Clear();
 140     delete fOutputList;
 141
 142   }
 143
 144   Printf("analysis done\n");
 145 }
 146
 147 //__________________________________________________________
 148 void AliHFEdisplacedElectrons::InitAnalysis(){
 149   //
 150   // init analysis (no intialization yet)
 151   //
 152
 153
 154   Printf("initialize analysis\n");
 155
 156 }
 157
 158
 159 //__________________________________________________________
 160 void AliHFEdisplacedElectrons::CreateOutputs(TList* const displacedList){
 161
 162   //
 163   //  create output fOutputList
 164   //
 165
 166   // THnSparseF
 167   // 8 interested electron sources: others, photon conv, direct photon,  pi0, eta, b, b->c, c
 168   // 21 possible DCA cuts XY
 169   // 21 possible DCA cuts Z
 170   // 13 pT bins
 171   // 10 rapidity bins
 172   // 10 azimuthal angle phi bins
 173
 174
 175   if(!displacedList) return;
 176
 177   fOutputList = displacedList;
 178   fOutputList -> SetName("information");
 179
 180
 181   // electron source
 182   Int_t nBinsEleSource = 8;
 183   Double_t minEleSource = -1.5;
 184   Double_t maxEleSource = 6.5;
 185   Double_t *binLimEleSource = new Double_t[nBinsEleSource+1];
 186   for(Int_t i=0; i<=nBinsEleSource; i++)
 187     binLimEleSource[i] = minEleSource + i*(maxEleSource-minEleSource)/nBinsEleSource;
 188
 189   // dca bins: the same as XY and Z
 190   // these should be variable bins as well
 191   Int_t nBinsDca = kNDcaMin-1;  // 12 bins
 192   Double_t minDca = -0.5;
 193   Double_t maxDca = 20.5;
 194   Double_t dcaBinWidth = (maxDca-minDca)/nBinsDca;
 195   Double_t *binLimDca = new Double_t[nBinsDca+1];
 196   for(Int_t i=0; i<=nBinsDca; i++)
 197     binLimDca[i] = minDca + i*dcaBinWidth;
 198
 199
 200  // pt bins
 201   Int_t nBinsPt = kNPtIntv-1;
 202   Double_t *binLimPt = new Double_t[nBinsPt+1];
 203   for(Int_t i=0; i<=nBinsPt; i++)
 204     binLimPt[i] = fgkPtIntv[i];  // variable bins
 205
 206   // rapidity bins
 207   Int_t nBinsRap = 10;
 208   Double_t minRap = -1.0;
 209   Double_t maxRap = 1.0;
 210   Double_t *binLimRap = new Double_t[nBinsRap+1];
 211   for(Int_t i=0; i<=nBinsRap; i++)
 212     binLimRap[i] = minRap + i*(maxRap-minRap)/nBinsRap;
 213
 214   // azumuthal phi angle
 215   Int_t nBinsPhi = 10;
 216   Double_t minPhi = 0;
 217   Double_t maxPhi = 2*TMath::Pi();
 218   Double_t *binLimPhi = new Double_t[nBinsPhi+1];
 219   for(Int_t i=0; i<=nBinsPhi; i++)
 220     binLimPhi[i] = minPhi + i*(maxPhi-minPhi)/nBinsPhi;
 221
 222
 223
 224   //for MC pionss
 225   const Int_t nVarPion = 5;
 226   Int_t iBinPion[nVarPion] = {nBinsDca, nBinsDca, nBinsPt, nBinsRap, nBinsPhi};
 227
 228   THnSparseF *fTHnSparseDcaMcPionInfo = NULL; // empty for the moment
 229   if(HasMCData()){
 230     fTHnSparseDcaMcPionInfo = new THnSparseF("dcaMcPionInfo",
 231                                              "MC info:;dcaXY [50 #mum];dcaZ [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 232                                              nVarPion, iBinPion);
 233     fTHnSparseDcaMcPionInfo->SetBinEdges(0, binLimDca);  // dca xy cut
 234     fTHnSparseDcaMcPionInfo->SetBinEdges(1, binLimDca);  // dca z cut
 235     fTHnSparseDcaMcPionInfo->SetBinEdges(2, binLimPt);   // pt
 236     fTHnSparseDcaMcPionInfo->SetBinEdges(3, binLimRap);  // rapidity
 237     fTHnSparseDcaMcPionInfo->SetBinEdges(4, binLimPhi);  // phi
 238     fTHnSparseDcaMcPionInfo->Sumw2();
 239
 240     fOutputList -> AddAt(fTHnSparseDcaMcPionInfo, kMCpion);
 241   }
 242
 243   // for MC electrons
 244   const Int_t nVar = 6;
 245   Int_t iBin[nVar] = {nBinsEleSource,nBinsDca, nBinsDca, nBinsPt, nBinsRap, nBinsPhi};
 246
 247   THnSparseF *fTHnSparseDcaMcEleInfo = NULL; // empty for the moment
 248   if(HasMCData()){
 249     fTHnSparseDcaMcEleInfo = new THnSparseF("dcaMcElectronInfo",
 250                                          "MC info:;ID [electron source id];dcaXY [50 #mum];dcaZ [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 251                                          nVar, iBin);
 252     fTHnSparseDcaMcEleInfo->SetBinEdges(0, binLimEleSource); // electron source
 253     fTHnSparseDcaMcEleInfo->SetBinEdges(1, binLimDca);  // dca xy cut
 254     fTHnSparseDcaMcEleInfo->SetBinEdges(2, binLimDca);  // dca z cut
 255     fTHnSparseDcaMcEleInfo->SetBinEdges(3, binLimPt);   // pt
 256     fTHnSparseDcaMcEleInfo->SetBinEdges(4, binLimRap);  // rapidity
 257     fTHnSparseDcaMcEleInfo->SetBinEdges(5, binLimPhi);  // phi
 258     fTHnSparseDcaMcEleInfo->Sumw2();
 259     fOutputList -> AddAt(fTHnSparseDcaMcEleInfo, kMCelectron);
 260   }
 261
 262
 263   // for ESD: HFE pid
 264
 265   THnSparseF *fTHnSparseDcaDataEleInfo = NULL;  // empty for the moment
 266   const Int_t nVarData = 5;
 267   Int_t iBinData[nVarData] = {nBinsDca, nBinsDca, nBinsPt, nBinsRap, nBinsPhi};
 268
 269   fTHnSparseDcaDataEleInfo = new THnSparseF("dcaDataElectronInfo",
 270                                          "Data info:;dcaXY [50 #mum];dcaZ [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 271                                          nVarData, iBinData);
 272   fTHnSparseDcaDataEleInfo->SetBinEdges(0, binLimDca);  // dca xy cut
 273   fTHnSparseDcaDataEleInfo->SetBinEdges(1, binLimDca);  // dca z cut
 274   fTHnSparseDcaDataEleInfo->SetBinEdges(2, binLimPt);   // pt
 275   fTHnSparseDcaDataEleInfo->SetBinEdges(3, binLimRap);  // rapidity
 276   fTHnSparseDcaDataEleInfo->SetBinEdges(4, binLimPhi);  // phi
 277   fTHnSparseDcaDataEleInfo->Sumw2();
 278
 279   fOutputList -> AddAt(fTHnSparseDcaDataEleInfo, kData);
 280
 281   AliInfo("THnSparse histograms are created\n");
 282   fOutputList->Print();
 283
 284 }
 285
 286
 287 //__________________________________________________________
 288 void AliHFEdisplacedElectrons::FillMCOutput(AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack, AliStack * const stack)
 289 {
 290
 291   // fill output
 292   if(!esdTrack) return;
 293   AliESDtrack *track = esdTrack;
 294   Double_t pt   = track->Pt();
 295   Double_t eta = track->Eta();
 296   Double_t phi = track->Phi();
 297
 298   Int_t label = track->GetLabel();
 299   if(label<0 || label>stack->GetNtrack()) return;
 300
 301   TParticle *particle = stack->Particle(label);
 302   if(!particle) return;
 303
 304
 305   // obtain impact parameters in xy and y
 306   const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex();
 307
 308   Float_t magneticField = 5;  // initialized as 5kG
 309   magneticField = fESDEvent->GetMagneticField();  // in kG
 310
 311
 312   Double_t dz[2];   // error of dca in cm
 313   Double_t covardz[3];
 314   track->PropagateToDCA(primVtx,magneticField, 1000., dz, covardz);
 315
 316
 317   Double_t dcaXY = TMath::Abs(dz[0])*1.0e4;  // conv dca in cm to dca in micron
 318   Double_t dcaZ = TMath::Abs(dz[1])*1.0e4;
 319
 320   // do PID with MC
 321
 322   if(HasMCData() && TMath::Abs(GetMCpid(stack, label))==kPDGelectron){
 323
 324     Int_t eleLabel = label;
 325
 326     const Int_t nvarMC=6;
 327     Double_t var[nvarMC];
 328     var[0] = -1;
 329     Int_t sourcePdg = ElectronFromSource(stack, eleLabel);
 330
 331     if(sourcePdg==kPDGgamma){ // check direct photon or not  // fixme
 332       if(ElePhotonDirect(stack, eleLabel)!=-1)
 333         var[0] = kEleDirectPhotonConv;
 334       else
 335         var[0] = kElePhotonConv;
 336       if(fDebugLevel>=10)
 337         printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 338     }   // photon or direct photon -> e
 339
 340     if(sourcePdg==kPDGpi0){
 341       var[0] = kElePi0;
 342       if(fDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 343     }
 344
 345     if(sourcePdg==kPDGeta){
 346       var[0] = kEleEta;
 347       if(fDebugLevel>=10)
 348         printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 349     }   // from eta -> e
 350
 351     if(TMath::Abs(sourcePdg%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 352        TMath::Abs(sourcePdg)/1000==kPDGbeauty ||
 353        TMath::Abs(sourcePdg)/100==kPDGbeauty ||
 354        TMath::Abs(sourcePdg)==kPDGbeauty){
 355       var[0]=kEleB;
 356       if(fDebugLevel>=10)
 357         printf("beauty======electron %d is from %d with id %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 358     } // direct beauty  -> e
 359     if(TMath::Abs(sourcePdg%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 1004**
 360        TMath::Abs(sourcePdg)/1000==kPDGcharm ||
 361        TMath::Abs(sourcePdg)/100==kPDGcharm ||
 362        TMath::Abs(sourcePdg)==kPDGcharm){
 363       // two cases:
 364       //            electron from b->c->e
 365       //            electron from c->e
 366       if(ElectronFromCharm(stack, eleLabel)!=-1){
 367         var[0] = ElectronFromCharm(stack, eleLabel);
 368         if(fDebugLevel>=10)
 369           printf("charm----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 370       }
 371     }  // charm electrons: b->c->e or c->e
 372
 373     if(fDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 374
 375     if(dcaXY<1000) var[1] = Int_t(dcaXY)/50;   // larger than 1mm should go to the last bin
 376     else
 377       var[1] = 20;
 378     if(dcaZ<1000) var[2] = Int_t(dcaZ)/50;
 379     else
 380       var[2] = 20;
 381
 382     var[3] = pt; // pt
 383     var[4] = eta; // eta
 384     var[5] = phi; // phi
 385
 386     (dynamic_cast<THnSparseF *>(fOutputList->At(kMCelectron)))->Fill(var);
 387   }
 388
 389   if(HasMCData() && TMath::Abs(GetMCpid(stack, label))==kPDGpion){
 390
 391     const Int_t nvarPionMC=5;
 392     Double_t varPion[nvarPionMC];
 393     if(dcaXY<1000)
 394       varPion[0] = Int_t(dcaXY)/50; // dca xy
 395     else
 396       varPion[0] = 20;
 397     if(dcaZ<1000)
 398       varPion[1] = Int_t(dcaZ)/50; // dca Z
 399     else
 400       varPion[1] = 20;
 401
 402 //     varPion[0] = TMath::Abs(dcaXY); // dca xy
 403 //     varPion[1] = TMath::Abs(dcaZ); // dca z
 404     varPion[2] = pt; // pt
 405     varPion[3] = eta; //eta
 406     varPion[4] = phi; // phi
 407
 408     (dynamic_cast<THnSparseF *>(fOutputList->At(kMCpion)))->Fill(varPion);
 409   }
 410
 411 }
 412
 413
 414
 415 //__________________________________________________________
 416 void AliHFEdisplacedElectrons::FillESDOutput(AliESDEvent * const fESD, AliESDtrack* const esdTrack)
 417 {
 418
 419   // fill output
 420   if(!esdTrack) return;
 421   AliESDtrack *track = esdTrack;
 422
 423   Double_t pt   = track->Pt();
 424   Double_t eta = track->Eta();
 425   Double_t phi = track->Phi();
 426
 427   // obtain impact parameters in xy and y
 428   const AliESDVertex *primVtx = fESD->GetPrimaryVertex();
 429
 430   Float_t magneticField = 5;  // initialized as 5kG
 431   magneticField = fESD->GetMagneticField();  // in kG
 432   Double_t dz[2];   // error of dca in cm
 433   Double_t covardz[3];
 434   track->PropagateToDCA(primVtx,magneticField, 1000., dz, covardz);
 435   Double_t dcaXY = TMath::Abs(dz[0]*1.0e4);  // conv dca in cm to dca in micron
 436   Double_t dcaZ = TMath::Abs(dz[1]*1.0e4);
 437
 438
 439   const Int_t nvarData = 5;
 440   Double_t varData[nvarData];
 441   // fixme
 442   if(dcaXY<1000)
 443     varData[0] = Int_t(dcaXY)/50; // dca xy
 444   else
 445     varData[0] = 20;
 446   if(dcaZ<1000)
 447     varData[1] = Int_t(dcaZ)/50; // dca Z
 448   else
 449     varData[1] = 20;
 450
 451   varData[2] = pt; // pt
 452   varData[3] = eta; //eta
 453   varData[4] = phi; // phi
 454
 455   (dynamic_cast<THnSparseF *>(fOutputList->At(kData)))->Fill(varData);
 456
 457 }
 458
 459 //__________________________________________________________
 460 Int_t AliHFEdisplacedElectrons::ElectronFromSource(AliStack * const stack, Int_t label) const
 461 {
 462
 463   //
 464   //  return electron source label via electron label
 465   //
 466
 467   // kEleSource is supposed to be either photon conv, direct photon conv, pi0, eta, beauty --> 0, 1, 2, 3, 4
 468
 469   Int_t eleLabel = label;
 470
 471   if(eleLabel<0 || eleLabel>stack->GetNtrack()) return -1;
 472
 473   TParticle *particle = stack->Particle(eleLabel);
 474   if(!particle) return -1;
 475
 476   Int_t motherLabel = particle->GetFirstMother();
 477   if(motherLabel<0 || motherLabel>stack->GetNtrack()) return -1;
 478   TParticle *motherPart = stack->Particle(motherLabel);
 479   if(!motherPart) return -1;
 480
 481   Int_t pdgCode = TMath::Abs(motherPart->GetPdgCode());
 482
 483   if(pdgCode==kPDGelectron) {
 484     if(fDebugLevel>=10) printf("particle label: %d...(motherLabel=%d : motherPdg=%d)  grandmother's pdg code was returned...%d \n",
 485                                label, motherLabel, pdgCode, ElectronFromSource(stack, motherLabel));
 486     return ElectronFromSource(stack, motherLabel);
 487   }
 488
 489   return pdgCode;
 490
 491 }
 492
 493 //__________________________________________________________
 494 Int_t AliHFEdisplacedElectrons::ElectronFromCharm(AliStack * const stack, Int_t label) const
 495 {
 496   //
 497   //  separate electron: kEleC from c->eX, kEleBC from b->c->eX
 498   //
 499
 500   Int_t motherLabel = label;
 501   TParticle *motherParticle = stack->Particle(motherLabel);  // mother part
 502   if(!motherParticle) return -1;
 503   Int_t gMotherLabel = motherParticle->GetFirstMother();  // grand mother
 504   if(gMotherLabel<0 || gMotherLabel>stack->GetNtrack()) return -1;
 505
 506   TParticle *gMotherPart = stack->Particle(gMotherLabel);
 507   if(!gMotherPart) return -1;
 508
 509   Int_t pdgCode = gMotherPart->GetPdgCode();
 510   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 511      TMath::Abs(pdgCode)/1000==kPDGbeauty || TMath::Abs(pdgCode)/100==kPDGbeauty || TMath::Abs(pdgCode)==kPDGbeauty) {
 512     if(fDebugLevel>=10)  printf("this electron label %d is from mother %d, and finally from %d\n", label, ElectronFromSource(stack, label), pdgCode);
 513     return kEleBC;
 514   }  // for sure it is from BC
 515
 516   else
 517     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||  // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 518        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 519        TMath::Abs(pdgCode)/100==kPDGcharm ||
 520        TMath::Abs(pdgCode)==kPDGcharm){
 521
 522       if(CharmFromBeauty(stack, gMotherLabel)!=-1)
 523         return kEleBC;
 524       else
 525         return kEleC;
 526
 527     }
 528
 529     else
 530       return -1;
 531
 532 }
 533
 534 //__________________________________________________________
 535 Int_t AliHFEdisplacedElectrons::CharmFromBeauty(AliStack * const stack, Int_t hfLabel) const
 536 {
 537
 538   //
 539   //  check if charm meson/hadron is from beauty decay
 540   //  -1, not from beauty
 541   //  return the label of the mother of this charm hadron
 542   //
 543
 544   Int_t charmLabel = hfLabel;
 545
 546   //  AliStack *stack = fMC->Stack();
 547   if(charmLabel<0 || charmLabel>stack->GetNtrack()) return -1;
 548   TParticle *particle = stack->Particle(charmLabel);
 549   if(!particle) return -1;
 550
 551   Int_t motherCharmLabel = particle->GetFirstMother();
 552   if(motherCharmLabel<0 || motherCharmLabel>stack->GetNtrack()) return -1;
 553
 554   TParticle *motherCharmPart = stack->Particle(motherCharmLabel);
 555   if(!motherCharmPart) return -1;
 556
 557   Int_t pdgCode = motherCharmPart->GetPdgCode();
 558
 559   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||   // for special intermediate meson states: like 100553
 560      TMath::Abs(pdgCode)/1000==kPDGbeauty ||
 561      TMath::Abs(pdgCode)/100==kPDGbeauty ||
 562      TMath::Abs(pdgCode)==kPDGbeauty)
 563     return motherCharmLabel;
 564   else
 565     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 566        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 567        TMath::Abs(pdgCode)/100==kPDGcharm)
 568       return CharmFromBeauty(stack, motherCharmLabel);   // loop over to see if charm is from beauty -- if yes, return the label of mother of the charm
 569
 570     else
 571       return -1;
 572 }
 573
 574 //__________________________________________________________
 575 Int_t AliHFEdisplacedElectrons::ElePhotonDirect(AliStack * const stack, Int_t label) const
 576 {
 577   //
 578   // electron is from photon, and check if this photon is direct
 579   //
 580
 581   Int_t eleLabel = label;
 582
 583   TParticle *particle = stack->Particle(eleLabel);
 584   if(particle->GetFirstMother()<0 || particle->GetFirstMother()>stack->GetNtrack())
 585     return -1;
 586
 587   TParticle *photonPart = stack->Particle(particle->GetFirstMother());
 588   if(!photonPart)
 589     return -1;
 590   Int_t motherPhotonLabel = photonPart->GetFirstMother();
 591   if(motherPhotonLabel<0 || motherPhotonLabel>stack->GetNtrack())
 592     return -1;
 593
 594   TParticle *motherPhotonPart = stack->Particle(motherPhotonLabel);
 595   if(!motherPhotonPart)
 596     return -1;
 597
 598   Int_t pdgMotherPhoton = motherPhotonPart->GetPdgCode();
 599   if(TMath::Abs(pdgMotherPhoton)<=10 || TMath::Abs(pdgMotherPhoton)==21)
 600     return 1;
 601
 602   else
 603     return -1;
 604
 605
 606 }
 607
 608 //__________________________________________________________
 609 Int_t AliHFEdisplacedElectrons::GetMCpid(AliStack* const stack, Int_t partLabel) const
 610 {
 611
 612   //
 613   // Simply pdg code
 614   //
 615
 616   Int_t label = partLabel;
 617 //   AliStack* stack = fMC->Stack();
 618   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 619
 620   // MC Information
 621   TParticle * particle = stack->Particle(label);
 622   if(!particle) return -1;
 623   Int_t pdg = particle->GetPdgCode();
 624
 625   return pdg;
 626
 627 }
 628
 629 //__________________________________________________________
 630 Int_t AliHFEdisplacedElectrons::GetMotherLabel(AliStack *const stack, Int_t eleLabel) const
 631 {
 632   //
 633   // Simply label of mother
 634   //
 635
 636
 637   Int_t label = eleLabel;
 638   //  AliStack* stack = fMC->Stack();
 639   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 640
 641   // MC Information
 642   TParticle * particle = stack->Particle(label);
 643   if(!particle) return -1;
 644
 645   return particle->GetFirstMother();
 646
 647 }
 648
 649
 650 //__________________________________________________________
 651 Float_t AliHFEdisplacedElectrons::GetRapidity(TParticle *part) const
 652 {
 653   // return rapidity
 654
 655   Float_t rapidity = -999;
 656   if((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)
 657     rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 658
 659   return rapidity;
 660 }
 661
 662
 663 //__________________________________________________________
 664 Float_t AliHFEdisplacedElectrons::GetTrackRapidity(AliESDtrack * const track) const
 665 {
 666   // return rapidity of electron
 667
 668   Float_t px = track->Px();
 669   Float_t py = track->Py();
 670   Float_t pz = track->Pz();
 671
 672   // electron mass 0.00051099906 GeV/c2
 673   TParticlePDG* electron = TDatabasePDG::Instance()->GetParticle(kPDGelectron);
 674   Double_t mass = electron->Mass();
 675
 676   Float_t en = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);
 677
 678   Float_t rapidity = -999;
 679   if((en - pz)*(en + pz)>0)
 680     rapidity = 0.5*(TMath::Log((en - pz)*(en + pz)));
 681
 682   return rapidity;
 683 }