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 //   Carlo Bombonati <Carlo.Bombonati@cern.ch>
  23 //
  24
  25 #include "TMath.h"
  26 #include "TList.h"
  27 #include "AliLog.h"
  28
  29 #include <TParticle.h>
  30 #include <TDatabasePDG.h>
  31 #include "THnSparse.h"
  32
  33 #include "AliMCEvent.h"
  34 #include "AliMCVertex.h"
  35 #include "AliMCParticle.h"
  36 #include "AliStack.h"
  37
  38 #include "AliESDEvent.h"
  39 #include "AliESDtrack.h"
  40
  41 #include "AliKFParticle.h"
  42 #include "AliKFVertex.h"
  43
  44 #include "AliVertexerTracks.h"
  45
  46
  47 #include "AliHFEdisplacedElectrons.h"
  48
  49 ClassImp(AliHFEdisplacedElectrons)
  50
  51 //__________________________________________________________
  52 const Float_t AliHFEdisplacedElectrons::fgkDcaMinPtIntv[13] = {
  53   //
  54   // define DCA low limits for single electrons cut in each Pt bin
  55   // preliminary numbers
  56   // these numbers should be replaced by the best numbers determined by
  57   // cut efficiency study: signal/background  (not used right now)
  58   //
  59   450, // 0.0 - 0.5
  60   450, // 0.5 - 1.0
  61   450, // 1.0 - 1.5
  62   500, // 1.5 - 2.0
  63   400, // 2.0 - 2.5
  64   300, // 2.5 - 3.0
  65   300, // 3.0 - 4.0
  66   300, // 4.0 - 5.0
  67   200, // 5.0 - 7.0
  68   150, // 7.0 - 9.0
  69   150, // 9.0 - 12.0
  70   100, //12.0 - 16.0
  71   50}; //16.0 - 20.0
  72 //__________________________________________________________
  73 const Float_t AliHFEdisplacedElectrons::fgkPtIntv[14] = {
  74   //
  75   // define pT bins for spectra of single electrons
  76   //
  77   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};
  78
  79 //__________________________________________________________
  80 const Char_t *AliHFEdisplacedElectrons::fgkKineVar[3] = {
  81   "y", "phi", "pt"};
  82
  83 //__________________________________________________________
  84 const Char_t *AliHFEdisplacedElectrons::fgkKineVarTitle[3] ={
  85   "rapidity;y;dN/dy",   "azimuthal;#phi;dN/d#phi",   "transverse momentum;p_{T};dN/dp_{T}",
  86 };
  87
  88 //__________________________________________________________
  89 AliHFEdisplacedElectrons::AliHFEdisplacedElectrons():
  90   fDeDebugLevel(0)
  91   , fNclustersITS(0)
  92   , fMinNprimVtxContributor(1)
  93   , fTHnSparseDcaMcEleInfo(NULL)
  94   , fTHnSparseDcaEsdEleInfo(NULL)
  95   , fTHnSparseDcaDataEleInfo(NULL)
  96   , fDeOutputList(0x0)
  97 {
  98   //
  99   // default constructor
 100   //
 101
 102 }
 103
 104 //__________________________________________________________
 105 AliHFEdisplacedElectrons::AliHFEdisplacedElectrons(const AliHFEdisplacedElectrons &ref):
 106   TObject(ref)
 107   , fDeDebugLevel(ref.fDeDebugLevel)
 108   , fNclustersITS(ref.fNclustersITS)
 109   , fMinNprimVtxContributor(ref.fMinNprimVtxContributor)
 110   , fTHnSparseDcaMcEleInfo(ref.fTHnSparseDcaMcEleInfo)
 111   , fTHnSparseDcaEsdEleInfo(ref.fTHnSparseDcaEsdEleInfo)
 112   , fTHnSparseDcaDataEleInfo(ref.fTHnSparseDcaDataEleInfo)
 113   , fDeOutputList(ref.fDeOutputList)
 114
 115 {
 116   //
 117   // copy constructor
 118   //
 119 }
 120
 121
 122 //__________________________________________________________
 123 AliHFEdisplacedElectrons&AliHFEdisplacedElectrons::operator=(const AliHFEdisplacedElectrons &ref)
 124 {
 125   //
 126   // Assignment operator
 127   //
 128
 129   if(this == &ref) return *this;
 130   AliHFEdisplacedElectrons::operator=(ref);
 131
 132   fDeDebugLevel = ref.fDeDebugLevel;
 133   fNclustersITS = ref.fNclustersITS;
 134   fMinNprimVtxContributor = ref.fMinNprimVtxContributor;
 135   fTHnSparseDcaMcEleInfo = ref.fTHnSparseDcaMcEleInfo;
 136   fTHnSparseDcaEsdEleInfo = ref.fTHnSparseDcaEsdEleInfo;
 137   fTHnSparseDcaDataEleInfo = ref.fTHnSparseDcaDataEleInfo;
 138   fDeOutputList = ref.fDeOutputList;
 139
 140   return *this;
 141 }
 142
 143 //__________________________________________________________
 144 AliHFEdisplacedElectrons::~AliHFEdisplacedElectrons()
 145 {
 146   //
 147   // default constructor
 148   //
 149
 150   if(fTHnSparseDcaMcEleInfo)
 151     delete fTHnSparseDcaMcEleInfo;
 152   if(fTHnSparseDcaEsdEleInfo)
 153     delete fTHnSparseDcaEsdEleInfo;
 154   if(fTHnSparseDcaDataEleInfo)
 155     delete fTHnSparseDcaDataEleInfo;
 156
 157   if(fDeOutputList){
 158     fDeOutputList->Clear();
 159     delete fDeOutputList;
 160
 161   }
 162
 163   Printf("analysis done\n");
 164 }
 165
 166 //__________________________________________________________
 167 void AliHFEdisplacedElectrons::InitAnalysis(){
 168   //
 169   // init analysis (no intialization yet)
 170   //
 171
 172
 173   Printf("initialize analysis\n");
 174
 175 }
 176
 177
 178 //__________________________________________________________
 179 void AliHFEdisplacedElectrons::CreateOutputs(TList* const displacedList){
 180
 181   //
 182   //  create output fDeOutputList
 183   //
 184
 185   // THnSparseF
 186   // 8+? interested electron sources: others, photon conv, direct photon,  pi0, eta, b, b->c, c + pion or missid, and missid pion
 187   // 41 possible DCA cuts XY:
 188   // 13 pT bins
 189   // 10 rapidity bins
 190   // 10 azimuthal angle phi bins
 191
 192
 193   if(!displacedList) return;
 194
 195   fDeOutputList = displacedList;
 196   fDeOutputList -> SetName("displacedElectrons");
 197
 198   // electron source
 199   Int_t nBinsEleSource = 10;
 200   Double_t minEleSource = -1.5;
 201   Double_t maxEleSource = 9.5;
 202   Double_t *binLimEleSource = new Double_t[nBinsEleSource+1];
 203   for(Int_t i=0; i<=nBinsEleSource; i++)
 204     binLimEleSource[i] = minEleSource + i*(maxEleSource-minEleSource)/nBinsEleSource;
 205
 206   // dca bins:  XY
 207   Int_t nBinsDcaXY = kNDcaMin-1;  // 41 bins
 208   Double_t minDcaXY = -20.5;
 209   Double_t maxDcaXY = 20.5;
 210   Double_t dcaXYBinWidth = (maxDcaXY-minDcaXY)/nBinsDcaXY;
 211   Double_t *binLimDcaXY = new Double_t[nBinsDcaXY+1];
 212   for(Int_t i=0; i<=nBinsDcaXY; i++)
 213     binLimDcaXY[i] = minDcaXY + i*dcaXYBinWidth;
 214
 215
 216   // dca bins:  Z
 217   Int_t nBinsDcaZ = kNDcaMin/2;  // 21 bins
 218   Double_t minDcaZ = -10.5;
 219   Double_t maxDcaZ = 10.5;
 220   Double_t dcaZBinWidth = (maxDcaZ-minDcaZ)/nBinsDcaZ;
 221   Double_t *binLimDcaZ = new Double_t[nBinsDcaZ+1];
 222   for(Int_t i=0; i<=nBinsDcaZ; i++)
 223     binLimDcaZ[i] = minDcaZ + i*dcaZBinWidth;
 224
 225  // pt bins
 226   Int_t nBinsPt = kNPtIntv-1;
 227   Double_t *binLimPt = new Double_t[nBinsPt+1];
 228   for(Int_t i=0; i<=nBinsPt; i++)
 229     binLimPt[i] = fgkPtIntv[i];  // variable bins
 230
 231   // rapidity bins
 232   Int_t nBinsRap = 10;
 233   Double_t minRap = -1.0;
 234   Double_t maxRap = 1.0;
 235   Double_t *binLimRap = new Double_t[nBinsRap+1];
 236   for(Int_t i=0; i<=nBinsRap; i++)
 237     binLimRap[i] = minRap + i*(maxRap-minRap)/nBinsRap;
 238
 239   // azumuthal phi angle
 240   Int_t nBinsPhi = 10;
 241   Double_t minPhi = 0;
 242   Double_t maxPhi = 2*TMath::Pi();
 243   Double_t *binLimPhi = new Double_t[nBinsPhi+1];
 244   for(Int_t i=0; i<=nBinsPhi; i++)
 245     binLimPhi[i] = minPhi + i*(maxPhi-minPhi)/nBinsPhi;
 246
 247
 248   fTHnSparseDcaMcEleInfo = 0x0;
 249   fTHnSparseDcaEsdEleInfo = 0x0;
 250   fTHnSparseDcaDataEleInfo = 0x0;
 251
 252   // for MC only: MC electron ID
 253   const Int_t nVarMc = 6;
 254   Int_t iBinMc[nVarMc] = {nBinsEleSource, nBinsDcaXY, nBinsDcaZ, nBinsPt, nBinsRap, nBinsPhi};
 255
 256   fTHnSparseDcaMcEleInfo = new THnSparseF("dcaMcElectronInfo",
 257                                            "MC electrons;electron source ID;mc dcaXY [50 #mum];mc dcaZ [100 #mum];mc pT [GeV/c];mc y [rapidity];mc #phi [rad];",
 258                                            nVarMc, iBinMc);
 259   fTHnSparseDcaMcEleInfo->SetBinEdges(0, binLimEleSource); // electron source
 260   fTHnSparseDcaMcEleInfo->SetBinEdges(1, binLimDcaXY);  // dca xy cut
 261   fTHnSparseDcaMcEleInfo->SetBinEdges(2, binLimDcaZ);  // dca z cut
 262   fTHnSparseDcaMcEleInfo->SetBinEdges(3, binLimPt);   // pt
 263   fTHnSparseDcaMcEleInfo->SetBinEdges(4, binLimRap);  // rapidity
 264   fTHnSparseDcaMcEleInfo->SetBinEdges(5, binLimPhi);  // phi
 265   fTHnSparseDcaMcEleInfo->Sumw2();
 266
 267   // for ESD with MC: HFE pid and MC pid
 268   const Int_t nVarEsd = 7;
 269   Int_t iBin[nVarEsd] = {nBinsEleSource,nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi};
 270
 271   fTHnSparseDcaEsdEleInfo
 272     = new THnSparseF("dcaEsdElectronInfo",
 273                      "ESD electrons;electron source ID;esd dcaXY [50 #mum];esd dcaZ [100 #mum];esd dcaXY KF [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 274                      nVarEsd, iBin);
 275
 276   fTHnSparseDcaEsdEleInfo->SetBinEdges(0, binLimEleSource); // electron source
 277   fTHnSparseDcaEsdEleInfo->SetBinEdges(1, binLimDcaXY);  // dca xy without current track
 278   fTHnSparseDcaEsdEleInfo->SetBinEdges(2, binLimDcaZ);  // dca z without current track
 279   fTHnSparseDcaEsdEleInfo->SetBinEdges(3, binLimDcaXY);  // dca xy kf without current track
 280   fTHnSparseDcaEsdEleInfo->SetBinEdges(4, binLimPt);   // pt esd
 281   fTHnSparseDcaEsdEleInfo->SetBinEdges(5, binLimRap);  // rapidity esd
 282   fTHnSparseDcaEsdEleInfo->SetBinEdges(6, binLimPhi);  // phi esd
 283   fTHnSparseDcaEsdEleInfo->Sumw2();
 284
 285   // for ESD data: HFE pid
 286   const Int_t nVarData = 6;
 287   Int_t iBinData[nVarData] = {nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi};
 288
 289   fTHnSparseDcaDataEleInfo
 290     = new THnSparseF("dcaDataElectronInfo",
 291                      "Data electrons;dcaXY [50 #mum];dcaZ [100 #mum];dcaXYkf [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 292                      nVarData, iBinData);
 293   fTHnSparseDcaDataEleInfo->SetBinEdges(0, binLimDcaXY);  // dca xy cut w/o
 294   fTHnSparseDcaDataEleInfo->SetBinEdges(1, binLimDcaZ);  // dca z cut w/o
 295   fTHnSparseDcaDataEleInfo->SetBinEdges(2, binLimDcaXY);  // dca xy kf cut
 296   fTHnSparseDcaDataEleInfo->SetBinEdges(3, binLimPt);   // pt
 297   fTHnSparseDcaDataEleInfo->SetBinEdges(4, binLimRap);  // rapidity
 298   fTHnSparseDcaDataEleInfo->SetBinEdges(5, binLimPhi);  // phi
 299   fTHnSparseDcaDataEleInfo->Sumw2();
 300
 301   fDeOutputList -> AddAt(fTHnSparseDcaMcEleInfo, kMcElectron);
 302   fDeOutputList -> AddAt(fTHnSparseDcaEsdEleInfo, kEsdElectron);
 303   fDeOutputList -> AddAt(fTHnSparseDcaDataEleInfo, kDataElectron);
 304
 305   AliInfo("THnSparse histograms are created\n");
 306   fDeOutputList->Print();
 307
 308 }
 309
 310
 311 //__________________________________________________________
 312 void AliHFEdisplacedElectrons::FillMcOutput(AliESDEvent *const fESD, AliMCEvent* const fMC, AliMCParticle* const mctrack)
 313 {
 314
 315   // fill output
 316   //0. after mc event cut
 317   //1. after checking stack, mcpart etc are valid
 318   //2. after PID
 319   //3. after event and track selection
 320
 321   AliStack *stack = fMC->Stack();
 322   TParticle *part = mctrack->Particle();
 323
 324   // obtain impact parameters in xy and z
 325   AliMCVertex *mcPrimVtx = (AliMCVertex *)fMC->GetPrimaryVertex();
 326   Double_t mcPrimV[3];
 327   mcPrimV[0] = mcPrimVtx->GetX();
 328   mcPrimV[1] = mcPrimVtx->GetY();
 329   mcPrimV[2] = mcPrimVtx->GetZ();
 330
 331  Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
 332
 333   Float_t vx = part->Vx();  // in cm
 334   Float_t vy = part->Vy();  // in cm
 335   Float_t vz = part->Vz();   // in cm
 336
 337   Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
 338
 339   Float_t mcpx = part->Px();
 340   Float_t mcpy = part->Py();
 341   Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
 342   Float_t mceta = part->Eta();
 343   Float_t mcphi = part->Phi();
 344
 345   Int_t pdg = part->GetPdgCode();
 346
 347   Int_t charge = 1;
 348   if(pdg==kPDGelectron || pdg==-kPDGpion) charge = -1;
 349
 350   // calculate mcDca ------------------------------------------------------------------
 351   const Float_t conv[2] = {1.783/1.6, 2.99792458};
 352   Float_t magneticField = 0;  // initialized as 5kG
 353   magneticField = fESD->GetMagneticField();  // in kG
 354   Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
 355
 356   Float_t radius;
 357   radius = TMath::Abs(radiusMc);
 358   Float_t nx = mcpx/mcpt;
 359   Float_t ny = mcpy/mcpt;
 360   Double_t dxy = vxy - mcVtxXY;   // in cm
 361   Double_t dvx = vx - mcPrimV[0]; // in cm
 362   Double_t dvy = vy - mcPrimV[1]; // in cm
 363   Float_t mcDcaXYm = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
 364   Double_t mcDca[2] = {mcDcaXYm*100, vz};  // in cm
 365   Double_t mcDcaXY = mcDca[0]*1.0e4;  // conv dca in cm to dca in micron
 366   Double_t mcDcaZ = mcDca[1]*1.0e4;
 367
 368   const Int_t nvarMC=6;
 369   Double_t var[nvarMC];
 370   var[0] = -1;
 371
 372   if(TMath::Abs(pdg)==kPDGelectron) {
 373
 374     Int_t eleLabel = mctrack->GetLabel();
 375     Int_t sourcePdg = ElectronFromSource(stack, eleLabel);
 376
 377     if(sourcePdg==kPDGgamma){ // check direct photon or not
 378       if(ElePhotonDirect(stack, eleLabel)!=-1)
 379         var[0] = kEleDirectPhotonConv;
 380       else
 381         var[0] = kElePhotonConv;
 382       if(fDeDebugLevel>=10)
 383         printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 384     }   // photon or direct photon -> e
 385
 386     if(sourcePdg==kPDGpi0){
 387       var[0] = kElePi0;
 388       if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 389     }
 390
 391     if(sourcePdg==kPDGeta){
 392       var[0] = kEleEta;
 393       if(fDeDebugLevel>=10)
 394         printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 395     }   // from eta -> e
 396
 397     if(TMath::Abs(sourcePdg%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 398        TMath::Abs(sourcePdg)/1000==kPDGbeauty ||
 399        TMath::Abs(sourcePdg)/100==kPDGbeauty ||
 400        TMath::Abs(sourcePdg)==kPDGbeauty){
 401       var[0]=kEleB;
 402       if(fDeDebugLevel>=10)
 403         printf("beauty======electron %d is from %d with id %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 404     } // direct beauty  -> e
 405     if(TMath::Abs(sourcePdg%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 1004**
 406        TMath::Abs(sourcePdg)/1000==kPDGcharm ||
 407        TMath::Abs(sourcePdg)/100==kPDGcharm ||
 408        TMath::Abs(sourcePdg)==kPDGcharm){
 409       // two cases:
 410       //            electron from b->c->e
 411       //            electron from c->e
 412       if(ElectronFromCharm(stack, eleLabel)!=-1){
 413         var[0] = ElectronFromCharm(stack, eleLabel);
 414         if(fDeDebugLevel>=10)
 415           printf("charm----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 416       }
 417     }  // charm electrons: b->c->e or c->e
 418
 419     if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 420   } // electron source endif
 421
 422   else
 423     if(TMath::Abs(pdg)==kPDGpion)
 424       var[0] = kPion;
 425
 426   if(TMath::Abs(mcDcaXY)<1000) var[1] = Int_t(mcDcaXY)/50;   // larger than 1mm should go to the last bin
 427   else
 428     var[1] = ((mcDcaXY>0)?1:-1)*20;
 429
 430   if(TMath::Abs(mcDcaZ)<1000) var[2] = Int_t(mcDcaZ)/100;
 431   else
 432     var[2] = ((mcDcaZ>0)?1:-1)*10;
 433
 434     var[3] = mcpt;  // pt
 435     var[4] = mceta; // eta
 436     var[5] = mcphi; // phi
 437
 438     (dynamic_cast<THnSparseF *>(fDeOutputList->At(kMcElectron)))->Fill(var);
 439 }
 440
 441
 442
 443 //__________________________________________________________
 444 void AliHFEdisplacedElectrons::FillEsdOutput(AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack, AliStack * const stack)
 445 {
 446   // after esd event selection, esd track cuts, and hfe pid
 447   // this is the case for ESD tracks, with MC information
 448
 449   AliESDtrack *track = esdTrack;
 450   Double_t pt   = track->Pt();
 451   Double_t eta = track->Eta();
 452   Double_t phi = track->Phi();
 453
 454   Int_t eleLabel = track->GetLabel();
 455   if(eleLabel<0 || eleLabel>stack->GetNtrack()) return;
 456
 457   TParticle *particle = stack->Particle(eleLabel);
 458   if(!particle) return;
 459
 460   // obtain impact parameters in xy and z
 461
 462   Float_t magneticField = 5;  // initialized as 5kG
 463   magneticField = fESDEvent->GetMagneticField();  // in kG
 464   Double_t beampiperadius=3.;
 465   const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex();
 466
 467
 468   const Int_t nvarESD = 7;
 469   Double_t var[nvarESD];
 470   var[0] = -1;
 471
 472   Int_t sourcePdg = -1;
 473
 474   if(TMath::Abs(particle->GetPdgCode())!=kPDGelectron)
 475     if(TMath::Abs(particle->GetPdgCode())!=kPDGpion)
 476       var[0] = kEleMissID;
 477     else
 478       var[0] = kEleMissIDpion;
 479
 480   else {
 481     sourcePdg = ElectronFromSource(stack, eleLabel);
 482
 483     if(sourcePdg==kPDGgamma){ // check direct photon or not
 484       if(ElePhotonDirect(stack, eleLabel)!=-1)
 485         var[0] = kEleDirectPhotonConv;
 486       else
 487         var[0] = kElePhotonConv;
 488       if(fDeDebugLevel>=10)
 489         printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 490     }   // photon or direct photon -> e
 491
 492     if(sourcePdg==kPDGpi0){
 493       var[0] = kElePi0;
 494       if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 495     }
 496
 497     if(sourcePdg==kPDGeta){
 498       var[0] = kEleEta;
 499       if(fDeDebugLevel>=10)
 500         printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 501     }   // from eta -> e
 502
 503     if(TMath::Abs(sourcePdg%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 504        TMath::Abs(sourcePdg)/1000==kPDGbeauty ||
 505        TMath::Abs(sourcePdg)/100==kPDGbeauty ||
 506        TMath::Abs(sourcePdg)==kPDGbeauty){
 507       var[0]=kEleB;
 508       if(fDeDebugLevel>=10)
 509         printf("beauty======electron %d is from %d with id %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 510     } // direct beauty  -> e
 511     if(TMath::Abs(sourcePdg%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 1004**
 512        TMath::Abs(sourcePdg)/1000==kPDGcharm ||
 513        TMath::Abs(sourcePdg)/100==kPDGcharm ||
 514        TMath::Abs(sourcePdg)==kPDGcharm){
 515       // two cases:
 516       //            electron from b->c->e
 517       //            electron from c->e
 518       if(ElectronFromCharm(stack, eleLabel)!=-1){
 519         var[0] = ElectronFromCharm(stack, eleLabel);
 520         if(fDeDebugLevel>=10)
 521           printf("charm----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 522       }
 523     }  // charm electrons: b->c->e or c->e
 524   }  // correct pid
 525
 526   // excluding current track
 527   // ---- beginning --- method from Andrea D 28.05.2010
 528   AliVertexerTracks *vertexer = new AliVertexerTracks(magneticField);
 529   vertexer->SetITSMode();
 530   vertexer->SetMinClusters(fNclustersITS);
 531   Int_t skipped[2];
 532   skipped[0] = (Int_t)track->GetID();
 533   vertexer->SetSkipTracks(1,skipped);
 534   AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent);
 535   delete vertexer; vertexer = NULL;
 536   if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return;
 537   // -- ending --- method from Andrea D 28.05.2010
 538
 539   Double_t dz[2];   // error of dca in cm
 540   Double_t covardz[3];
 541   if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
 542
 543   Double_t dcaXY = dz[0]*1.0e4;  // conv dca in cm to dca in micron
 544   Double_t dcaZ = dz[1]*1.0e4;
 545
 546   if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 547
 548   if(TMath::Abs(dcaXY)<1000) var[1] = Int_t(dcaXY)/50;   // larger than 1mm should go to the last bin
 549   else
 550     var[1] = ((dcaXY>0)?1:-1)*20;
 551
 552   if(TMath::Abs(dcaZ)<1000) var[2] = Int_t(dcaZ)/100;
 553   else
 554     var[2] = ((dcaZ>0)?1:-1)*10;
 555
 556   // calculate dca using AliKFParticle class------------------------------------------------------------------
 557   Float_t  kfDcaXY = 0;
 558   Int_t trkID = track->GetID();
 559   AliKFParticle::SetField(magneticField);
 560   AliKFParticle kfParticle(*track, particle->GetPdgCode());
 561   // prepare kfprimary vertex
 562   AliKFVertex kfESDprimary;
 563   // Reconstruct Primary Vertex (with ESD tracks)
 564   Int_t n=primVtx->GetNIndices();
 565   if (n>0 && primVtx->GetStatus()){
 566     kfESDprimary = AliKFVertex(*primVtx);
 567     UShort_t *priIndex = primVtx->GetIndices();
 568     for (Int_t i=0;i<n;i++){
 569       Int_t idx = Int_t(priIndex[i]);
 570       if (idx == trkID){
 571         kfESDprimary -= kfParticle;
 572         kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4;
 573       }  // remove current track from this calculation
 574     }  // loop over all primary vertex contributors
 575   }
 576   // end of KF dca
 577
 578   if(TMath::Abs(kfDcaXY)<1000)
 579     var[3] = Int_t(kfDcaXY)/50;
 580   else
 581     var[3] = (kfDcaXY>0?1:-1)*20;
 582
 583   var[4] = pt;  // pt
 584   var[5] = eta; // eta
 585   var[6] = phi; // phi
 586
 587   (dynamic_cast<THnSparseF *>(fDeOutputList->At(kEsdElectron)))->Fill(var);
 588
 589 }
 590
 591 //__________________________________________________________
 592 void AliHFEdisplacedElectrons::FillDataOutput(AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack)
 593 {
 594
 595   // this is pure data, without MC information at all
 596   // fill output: with HFE pid selection of electrons after all track quality cuts
 597
 598   if(!esdTrack) return;
 599   AliESDtrack *track = esdTrack;
 600
 601   Double_t pt   = track->Pt();
 602   Double_t eta = track->Eta();
 603   Double_t phi = track->Phi();
 604
 605   // obtain impact parameters in xy and y
 606   const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex();
 607
 608   Float_t magneticField = 5;  // initialized as 5kG
 609   magneticField = fESDEvent->GetMagneticField();  // in kG
 610   Double_t beampiperadius=3.;
 611
 612
 613   const Int_t nvarData = 6;
 614   Double_t varData[nvarData];
 615
 616   //
 617   // excluding current track
 618   //
 619
 620   //------ beginning --- method from Andrea D 28.05.2010
 621   AliVertexerTracks *vertexer = new AliVertexerTracks(fESDEvent->GetMagneticField());
 622   vertexer->SetITSMode();
 623   vertexer->SetMinClusters(fNclustersITS);
 624   Int_t skipped[2];
 625   skipped[0] = (Int_t)track->GetID();
 626   vertexer->SetSkipTracks(1,skipped);
 627   AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent);
 628   delete vertexer; vertexer = NULL;
 629   if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return;
 630   //------ ending --- method from Andrea D 28.05.2010
 631
 632   Double_t dz[2];   // error of dca in cm
 633   Double_t covardz[3];
 634
 635   if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
 636
 637   Double_t dcaXY = dz[0]*1.0e4;  // conv dca in cm to dca in micron
 638   Double_t dcaZ = dz[1]*1.0e4;
 639
 640   if(TMath::Abs(dcaXY)<1000) varData[0] = Int_t(dcaXY)/50;   // larger than 1mm should go to the last bin
 641   else
 642     varData[0] = (dcaXY>0?1:-1)*20;
 643   if(TMath::Abs(dcaZ)<1000) varData[1] = Int_t(dcaZ)/100;
 644   else
 645     varData[1] = (dcaZ>0?1:-1)*10;
 646
 647
 648   // calculate dca using AliKFParticle class------------------------------------------------------------------
 649   Float_t  kfDcaXY = 0;
 650   Int_t trkID = track->GetID();
 651   AliKFParticle::SetField(magneticField);
 652   AliKFParticle kfParticle(*track, -11*track->Charge());
 653   // prepare kfprimary vertex
 654   AliKFVertex kfESDprimary;
 655   // Reconstruct Primary Vertex (with ESD tracks)
 656   Int_t n=primVtx->GetNIndices();
 657   if (n>0 && primVtx->GetStatus()){
 658     kfESDprimary = AliKFVertex(*primVtx);
 659     UShort_t *priIndex = primVtx->GetIndices();
 660     for (Int_t i=0;i<n;i++){
 661       Int_t idx = Int_t(priIndex[i]);
 662       if (idx == trkID){
 663           kfESDprimary -= kfParticle;
 664           kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4;
 665       }  // remove current track from this calculation
 666     }  // loop over all primary vertex contributors
 667   }
 668   if(TMath::Abs(kfDcaXY)<1000)
 669     varData[2] = Int_t(kfDcaXY)/50;
 670   else
 671     varData[2] = ((kfDcaXY)>0?1:-1)*20;
 672   // end of KF dca
 673
 674   varData[3] = pt; // pt
 675   varData[4] = eta; //eta
 676   varData[5] = phi; // phi
 677
 678   (dynamic_cast<THnSparseF *>(fDeOutputList->At(kDataElectron)))->Fill(varData);
 679
 680 }
 681
 682 //__________________________________________________________
 683 Int_t AliHFEdisplacedElectrons::ElectronFromSource(AliStack * const stack, Int_t label) const
 684 {
 685
 686   //
 687   //  return electron source label via electron label
 688   //
 689
 690   // kEleSource is supposed to be either photon conv, direct photon conv, pi0, eta, beauty --> 0, 1, 2, 3, 4
 691
 692   Int_t eleLabel = label;
 693
 694   if(eleLabel<0 || eleLabel>stack->GetNtrack()) return -1;
 695
 696   TParticle *particle = stack->Particle(eleLabel);
 697   if(!particle) return -1;
 698
 699   Int_t motherLabel = particle->GetFirstMother();
 700   if(motherLabel<0 || motherLabel>stack->GetNtrack()) return -1;
 701   TParticle *motherPart = stack->Particle(motherLabel);
 702   if(!motherPart) return -1;
 703
 704   Int_t pdgCode = TMath::Abs(motherPart->GetPdgCode());
 705
 706   if(pdgCode==kPDGelectron) {
 707     if(fDeDebugLevel>=10) printf("particle label: %d...(motherLabel=%d : motherPdg=%d)  grandmother's pdg code was returned...%d \n",
 708                                label, motherLabel, pdgCode, ElectronFromSource(stack, motherLabel));
 709     return ElectronFromSource(stack, motherLabel);
 710   }
 711
 712   return pdgCode;
 713
 714 }
 715
 716 //__________________________________________________________
 717 Int_t AliHFEdisplacedElectrons::ElectronFromCharm(AliStack * const stack, Int_t label) const
 718 {
 719   //
 720   //  separate electron: kEleC from c->eX, kEleBC from b->c->eX
 721   //
 722
 723   Int_t motherLabel = label;
 724   TParticle *motherParticle = stack->Particle(motherLabel);  // mother part
 725   if(!motherParticle) return -1;
 726   Int_t gMotherLabel = motherParticle->GetFirstMother();  // grand mother
 727   if(gMotherLabel<0 || gMotherLabel>stack->GetNtrack()) return -1;
 728
 729   TParticle *gMotherPart = stack->Particle(gMotherLabel);
 730   if(!gMotherPart) return -1;
 731
 732   Int_t pdgCode = gMotherPart->GetPdgCode();
 733   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 734      TMath::Abs(pdgCode)/1000==kPDGbeauty || TMath::Abs(pdgCode)/100==kPDGbeauty || TMath::Abs(pdgCode)==kPDGbeauty) {
 735     if(fDeDebugLevel>=10)  printf("this electron label %d is from mother %d, and finally from %d\n", label, ElectronFromSource(stack, label), pdgCode);
 736     return kEleBC;
 737   }  // for sure it is from BC
 738
 739   else
 740     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||  // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 741        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 742        TMath::Abs(pdgCode)/100==kPDGcharm ||
 743        TMath::Abs(pdgCode)==kPDGcharm){
 744
 745       if(CharmFromBeauty(stack, gMotherLabel)!=-1)
 746         return kEleBC;
 747       else
 748         return kEleC;
 749
 750     }
 751
 752     else
 753       return -1;
 754
 755 }
 756
 757 //__________________________________________________________
 758 Int_t AliHFEdisplacedElectrons::CharmFromBeauty(AliStack * const stack, Int_t hfLabel) const
 759 {
 760
 761   //
 762   //  check if charm meson/hadron is from beauty decay
 763   //  -1, not from beauty
 764   //  return the label of the mother of this charm hadron
 765   //
 766
 767   Int_t charmLabel = hfLabel;
 768
 769   //  AliStack *stack = fMC->Stack();
 770   if(charmLabel<0 || charmLabel>stack->GetNtrack()) return -1;
 771   TParticle *particle = stack->Particle(charmLabel);
 772   if(!particle) return -1;
 773
 774   Int_t motherCharmLabel = particle->GetFirstMother();
 775   if(motherCharmLabel<0 || motherCharmLabel>stack->GetNtrack()) return -1;
 776
 777   TParticle *motherCharmPart = stack->Particle(motherCharmLabel);
 778   if(!motherCharmPart) return -1;
 779
 780   Int_t pdgCode = motherCharmPart->GetPdgCode();
 781
 782   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||   // for special intermediate meson states: like 100553
 783      TMath::Abs(pdgCode)/1000==kPDGbeauty ||
 784      TMath::Abs(pdgCode)/100==kPDGbeauty ||
 785      TMath::Abs(pdgCode)==kPDGbeauty)
 786     return motherCharmLabel;
 787   else
 788     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 789        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 790        TMath::Abs(pdgCode)/100==kPDGcharm)
 791       return CharmFromBeauty(stack, motherCharmLabel);   // loop over to see if charm is from beauty -- if yes, return the label of mother of the charm
 792
 793     else
 794       return -1;
 795 }
 796
 797 //__________________________________________________________
 798 Int_t AliHFEdisplacedElectrons::ElePhotonDirect(AliStack * const stack, Int_t label) const
 799 {
 800   //
 801   // electron is from photon, and check if this photon is direct
 802   //
 803
 804   Int_t eleLabel = label;
 805
 806   TParticle *particle = stack->Particle(eleLabel);
 807   if(particle->GetFirstMother()<0 || particle->GetFirstMother()>stack->GetNtrack())
 808     return -1;
 809
 810   TParticle *photonPart = stack->Particle(particle->GetFirstMother());
 811   if(!photonPart)
 812     return -1;
 813   Int_t motherPhotonLabel = photonPart->GetFirstMother();
 814   if(motherPhotonLabel<0 || motherPhotonLabel>stack->GetNtrack())
 815     return -1;
 816
 817   TParticle *motherPhotonPart = stack->Particle(motherPhotonLabel);
 818   if(!motherPhotonPart)
 819     return -1;
 820
 821   Int_t pdgMotherPhoton = motherPhotonPart->GetPdgCode();
 822   if(TMath::Abs(pdgMotherPhoton)<=10 || TMath::Abs(pdgMotherPhoton)==21)
 823     return 1;
 824
 825   else
 826     return -1;
 827
 828
 829 }
 830
 831 //__________________________________________________________
 832 Int_t AliHFEdisplacedElectrons::GetMCpid(AliStack* const stack, Int_t partLabel) const
 833 {
 834
 835   //
 836   // Simply pdg code
 837   //
 838
 839   Int_t label = partLabel;
 840 //   AliStack* stack = fMC->Stack();
 841   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 842
 843   // MC Information
 844   TParticle * particle = stack->Particle(label);
 845   if(!particle) return -1;
 846   Int_t pdg = particle->GetPdgCode();
 847
 848   return pdg;
 849
 850 }
 851
 852 //__________________________________________________________
 853 Int_t AliHFEdisplacedElectrons::GetMotherLabel(AliStack *const stack, Int_t eleLabel) const
 854 {
 855   //
 856   // Simply label of mother
 857   //
 858
 859
 860   Int_t label = eleLabel;
 861   //  AliStack* stack = fMC->Stack();
 862   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 863
 864   // MC Information
 865   TParticle * particle = stack->Particle(label);
 866   if(!particle) return -1;
 867
 868   return particle->GetFirstMother();
 869
 870 }
 871
 872
 873 //__________________________________________________________
 874 Float_t AliHFEdisplacedElectrons::GetRapidity(TParticle *part) const
 875 {
 876   // return rapidity
 877
 878   Float_t rapidity = -999;
 879   if((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)
 880     rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 881
 882   return rapidity;
 883 }
 884
 885
 886 //__________________________________________________________
 887 Float_t AliHFEdisplacedElectrons::GetTrackRapidity(AliESDtrack * const track) const
 888 {
 889   // return rapidity of electron
 890
 891   Float_t px = track->Px();
 892   Float_t py = track->Py();
 893   Float_t pz = track->Pz();
 894
 895   // electron mass 0.00051099906 GeV/c2
 896   TParticlePDG* electron = TDatabasePDG::Instance()->GetParticle(kPDGelectron);
 897   Double_t mass = electron->Mass();
 898
 899   Float_t en = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);
 900
 901   Float_t rapidity = -999;
 902   if((en - pz)*(en + pz)>0)
 903     rapidity = 0.5*(TMath::Log((en - pz)*(en + pz)));
 904
 905   return rapidity;
 906 }