PWGHF/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 = 8.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   // production radii
 248   Int_t nBinsR = 30;
 249   Double_t minR = 0;
 250   Double_t maxR = 15;
 251   Double_t *binLimR = new Double_t[nBinsR+1];
 252   for(Int_t i=0; i<=nBinsR; i++)
 253     binLimR[i] = minR + i*(maxR-minR)/nBinsR;
 254
 255   // status
 256   Int_t nBinsStat = 11;
 257   Double_t minStat = 0;
 258   Double_t maxStat = 11;
 259   Double_t *binLimStat = new Double_t[nBinsStat+1];
 260   for(Int_t i=0; i<=nBinsStat; i++)
 261     binLimStat[i] = minStat + i*(maxStat-minStat)/nBinsStat;
 262
 263   fTHnSparseDcaMcEleInfo = 0x0;
 264   fTHnSparseDcaEsdEleInfo = 0x0;
 265   fTHnSparseDcaDataEleInfo = 0x0;
 266
 267   // for MC only: MC electron ID
 268   const Int_t nVarMc = 8;
 269   Int_t iBinMc[nVarMc] = {nBinsEleSource, nBinsDcaXY, nBinsDcaZ, nBinsPt, nBinsRap, nBinsPhi, nBinsR, nBinsStat};
 270
 271   fTHnSparseDcaMcEleInfo = new THnSparseF("dcaMcElectronInfo",
 272                                            "MC electrons;electron source ID;mc dcaXY [50 #mum];mc dcaZ [100 #mum];mc pT [GeV/c];mc y [rapidity];mc #phi [rad];mc R [cm];Status Code",
 273                                            nVarMc, iBinMc);
 274
 275   fTHnSparseDcaMcEleInfo->SetBinEdges(0, binLimEleSource); // electron source
 276   fTHnSparseDcaMcEleInfo->SetBinEdges(1, binLimDcaXY);  // dca xy cut
 277   fTHnSparseDcaMcEleInfo->SetBinEdges(2, binLimDcaZ);  // dca z cut
 278   fTHnSparseDcaMcEleInfo->SetBinEdges(3, binLimPt);   // pt
 279   fTHnSparseDcaMcEleInfo->SetBinEdges(4, binLimRap);  // rapidity
 280   fTHnSparseDcaMcEleInfo->SetBinEdges(5, binLimPhi);  // phi
 281   fTHnSparseDcaMcEleInfo->SetBinEdges(6, binLimR);
 282   fTHnSparseDcaMcEleInfo->SetBinEdges(7, binLimStat);
 283   fTHnSparseDcaMcEleInfo->Sumw2();
 284
 285   // for ESD with MC: HFE pid and MC pid
 286   const Int_t nVarEsd = 9;
 287   Int_t iBin[nVarEsd] = {nBinsEleSource,nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi, nBinsR, nBinsStat};
 288
 289   fTHnSparseDcaEsdEleInfo
 290     = new THnSparseF("dcaEsdElectronInfo",
 291                      "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];R [cm];Status Code",
 292                      nVarEsd, iBin);
 293
 294   fTHnSparseDcaEsdEleInfo->SetBinEdges(0, binLimEleSource); // electron source
 295   fTHnSparseDcaEsdEleInfo->SetBinEdges(1, binLimDcaXY);  // dca xy without current track
 296   fTHnSparseDcaEsdEleInfo->SetBinEdges(2, binLimDcaZ);  // dca z without current track
 297   fTHnSparseDcaEsdEleInfo->SetBinEdges(3, binLimDcaXY);  // dca xy kf without current track
 298   fTHnSparseDcaEsdEleInfo->SetBinEdges(4, binLimPt);   // pt esd
 299   fTHnSparseDcaEsdEleInfo->SetBinEdges(5, binLimRap);  // rapidity esd
 300   fTHnSparseDcaEsdEleInfo->SetBinEdges(6, binLimPhi);  // phi esd
 301   fTHnSparseDcaEsdEleInfo->SetBinEdges(7, binLimR);
 302   fTHnSparseDcaEsdEleInfo->SetBinEdges(8, binLimStat);
 303   fTHnSparseDcaEsdEleInfo->Sumw2();
 304
 305
 306   // for ESD data: HFE pid
 307   const Int_t nVarData = 6;
 308   Int_t iBinData[nVarData] = {nBinsDcaXY, nBinsDcaZ, nBinsDcaXY, nBinsPt, nBinsRap, nBinsPhi};
 309
 310   fTHnSparseDcaDataEleInfo
 311     = new THnSparseF("dcaDataElectronInfo",
 312                      "Data electrons;dcaXY [50 #mum];dcaZ [100 #mum];dcaXYkf [50 #mum];pT [GeV/c];y [rapidity];#phi [rad];",
 313                      nVarData, iBinData);
 314   fTHnSparseDcaDataEleInfo->SetBinEdges(0, binLimDcaXY);  // dca xy cut w/o
 315   fTHnSparseDcaDataEleInfo->SetBinEdges(1, binLimDcaZ);  // dca z cut w/o
 316   fTHnSparseDcaDataEleInfo->SetBinEdges(2, binLimDcaXY);  // dca xy kf cut
 317   fTHnSparseDcaDataEleInfo->SetBinEdges(3, binLimPt);   // pt
 318   fTHnSparseDcaDataEleInfo->SetBinEdges(4, binLimRap);  // rapidity
 319   fTHnSparseDcaDataEleInfo->SetBinEdges(5, binLimPhi);  // phi
 320   fTHnSparseDcaDataEleInfo->Sumw2();
 321
 322   fDeOutputList -> AddAt(fTHnSparseDcaMcEleInfo, kMcElectron);
 323   fDeOutputList -> AddAt(fTHnSparseDcaEsdEleInfo, kEsdElectron);
 324   fDeOutputList -> AddAt(fTHnSparseDcaDataEleInfo, kDataElectron);
 325
 326   AliInfo("THnSparse histograms are created\n");
 327   fDeOutputList->Print();
 328
 329 }
 330
 331
 332 //__________________________________________________________
 333 void AliHFEdisplacedElectrons::FillMcOutput(const AliESDEvent *const fESD, AliMCEvent* const fMC, const AliMCParticle* const mctrack)
 334 {
 335
 336   // fill output
 337   //0. after mc event cut
 338   //1. after checking stack, mcpart etc are valid
 339   //2. after PID
 340   //3. after event and track selection
 341
 342   AliStack *stack = fMC->Stack();
 343   TParticle *part = mctrack->Particle();
 344
 345   // obtain impact parameters in xy and z
 346   AliMCVertex *mcPrimVtx = (AliMCVertex *)fMC->GetPrimaryVertex();
 347   Double_t mcPrimV[3];
 348   mcPrimV[0] = mcPrimVtx->GetX();
 349   mcPrimV[1] = mcPrimVtx->GetY();
 350   mcPrimV[2] = mcPrimVtx->GetZ();
 351
 352   Double_t mcVtxXY = TMath::Abs(mcPrimV[0]*mcPrimV[0] + mcPrimV[1]*mcPrimV[1]);
 353
 354   Float_t vx = part->Vx();  // in cm
 355   Float_t vy = part->Vy();  // in cm
 356   Float_t vz = part->Vz();   // in cm
 357
 358   Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
 359
 360   Float_t mcpx = part->Px();
 361   Float_t mcpy = part->Py();
 362   Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
 363   Float_t mceta = part->Eta();
 364   Float_t mcphi = part->Phi();
 365
 366   // begin new stuff
 367   Double_t mcR = part->R();
 368   Int_t mcStatus = part->GetStatusCode();
 369   // end new staff
 370
 371   Int_t pdg = part->GetPdgCode();
 372
 373   Int_t charge = 1;
 374   if(pdg==kPDGelectron || pdg==-kPDGpion) charge = -1;
 375
 376   // calculate mcDca ------------------------------------------------------------------
 377   const Float_t conv[2] = {1.783/1.6, 2.99792458};
 378   Float_t magneticField = 0;  // initialized as 5kG
 379   magneticField = fESD->GetMagneticField();  // in kG
 380   Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla, radius in meter
 381
 382   Float_t radius;
 383   radius = TMath::Abs(radiusMc);
 384   Float_t nx = mcpx/mcpt;
 385   Float_t ny = mcpy/mcpt;
 386   Double_t dxy = vxy - mcVtxXY;   // in cm
 387   Double_t dvx = vx - mcPrimV[0]; // in cm
 388   Double_t dvy = vy - mcPrimV[1]; // in cm
 389   Float_t mcDcaXYm = (radius - TMath::Sqrt(dxy*dxy/100./100. + radius*radius + 2*radius*charge*(dvx*ny-dvy*nx)/100.)) ;  // in meters
 390   Double_t mcDca[2] = {mcDcaXYm*100, vz};  // in cm
 391   Double_t mcDcaXY = mcDca[0]*1.0e4;  // conv dca in cm to dca in micron
 392   Double_t mcDcaZ = mcDca[1]*1.0e4;
 393
 394   const Int_t nvarMC=8;
 395   Double_t var[nvarMC];
 396   var[0] = -1;
 397
 398   if(TMath::Abs(pdg)==kPDGelectron) {
 399
 400     Int_t eleLabel = mctrack->GetLabel();
 401     Int_t sourcePdg = ElectronFromSource(stack, eleLabel);
 402
 403     if(sourcePdg==kPDGgamma){ // check direct photon or not
 404       if(ElePhotonDirect(stack, eleLabel)!=-1)
 405         var[0] = kEleDirectPhotonConv;
 406       else
 407         var[0] = kElePhotonConv;
 408       if(fDeDebugLevel>=10)
 409         printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 410     }   // photon or direct photon -> e
 411
 412     if(sourcePdg==kPDGpi0){
 413       var[0] = kElePi0;
 414       if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 415     }
 416
 417     if(sourcePdg==kPDGeta){
 418       var[0] = kEleEta;
 419       if(fDeDebugLevel>=10)
 420         printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 421     }   // from eta -> e
 422
 423     if(TMath::Abs(sourcePdg%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 424        TMath::Abs(sourcePdg)/1000==kPDGbeauty ||
 425        TMath::Abs(sourcePdg)/100==kPDGbeauty ||
 426        TMath::Abs(sourcePdg)==kPDGbeauty){
 427       var[0]=kEleB;
 428       if(fDeDebugLevel>=10)
 429         printf("beauty======electron %d is from %d with id %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 430     } // direct beauty  -> e
 431     if(TMath::Abs(sourcePdg%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 1004**
 432        TMath::Abs(sourcePdg)/1000==kPDGcharm ||
 433        TMath::Abs(sourcePdg)/100==kPDGcharm ||
 434        TMath::Abs(sourcePdg)==kPDGcharm){
 435       // two cases:
 436       //            electron from b->c->e
 437       //            electron from c->e
 438       if(ElectronFromCharm(stack, eleLabel)!=-1){
 439         var[0] = ElectronFromCharm(stack, eleLabel);
 440         if(fDeDebugLevel>=10)
 441           printf("charm----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 442       }
 443     }  // charm electrons: b->c->e or c->e
 444
 445     if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 446   } // electron source endif
 447
 448   else
 449     if(TMath::Abs(pdg)==kPDGpion)
 450       var[0] = kPion;
 451
 452   if(TMath::Abs(mcDcaXY)<1000) var[1] = Int_t(mcDcaXY)/50;   // larger than 1mm should go to the last bin
 453   else
 454     var[1] = ((mcDcaXY>0)?1:-1)*20;
 455
 456   if(TMath::Abs(mcDcaZ)<1000) var[2] = Int_t(mcDcaZ)/100;
 457   else
 458     var[2] = ((mcDcaZ>0)?1:-1)*10;
 459
 460     var[3] = mcpt;      // pt
 461     var[4] = mceta;     // eta
 462     var[5] = mcphi;     // phi
 463     var[6] = mcR;       // production radius
 464     var[7] = mcStatus;  // internal status
 465
 466
 467     (static_cast<THnSparseF *>(fDeOutputList->At(kMcElectron)))->Fill(var);
 468 }
 469
 470
 471
 472 //__________________________________________________________
 473 void AliHFEdisplacedElectrons::FillEsdOutput(const AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack, AliStack * const stack)
 474 {
 475   // after esd event selection, esd track cuts, and hfe pid
 476   // this is the case for ESD tracks, with MC information
 477
 478   AliESDtrack *track = esdTrack;
 479   Double_t pt  = track->Pt();
 480   Double_t eta = track->Eta();
 481   Double_t phi = track->Phi();
 482
 483   Int_t eleLabel = track->GetLabel();
 484   if(eleLabel<0 || eleLabel>stack->GetNtrack()) return;
 485
 486   TParticle *particle = stack->Particle(eleLabel);
 487   if(!particle) return;
 488
 489   // begin new stuff
 490   Double_t mcR = particle->R();
 491   Int_t mcStatus = particle->GetStatusCode();
 492   // end new staff
 493
 494   // obtain impact parameters in xy and z
 495   Float_t magneticField = 5;  // initialized as 5kG
 496   magneticField = fESDEvent->GetMagneticField();  // in kG
 497   Double_t beampiperadius=3.;
 498   const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex();
 499
 500
 501   const Int_t nvarESD = 9;
 502   Double_t var[nvarESD];
 503   var[0] = -1;
 504
 505   Int_t sourcePdg = -1;
 506
 507   if(TMath::Abs(particle->GetPdgCode())==kPDGelectron){
 508
 509     sourcePdg = ElectronFromSource(stack, eleLabel);
 510
 511     if(sourcePdg==kPDGgamma){ // check direct photon or not
 512       if(ElePhotonDirect(stack, eleLabel)!=-1)
 513         var[0] = kEleDirectPhotonConv;
 514       else
 515         var[0] = kElePhotonConv;
 516       if(fDeDebugLevel>=10)
 517         printf("photonconv=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 518     } else if(sourcePdg==kPDGpi0){ // check dalitz
 519       var[0] = kElePi0;
 520       if(fDeDebugLevel>=10) printf("pi0======this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 521     } else if(sourcePdg==kPDGeta){ // check eta
 522       var[0] = kEleEta;
 523       if(fDeDebugLevel>=10)
 524         printf("eta=====this electron %d is from %d, source id=%.1f\n", eleLabel, sourcePdg, var[0]);
 525     } else if(IsB(sourcePdg)) var[0] = kEleB; // check beauty
 526       else if(IsC(sourcePdg)) var[0] = CheckCharm(stack,eleLabel); // check charm
 527
 528   } else if(TMath::Abs(particle->GetPdgCode())==kPDGpion) var[0] = kEleMissIDpion;
 529   else var[0] = kEleMissID;
 530   // ---- PID END ----
 531
 532
 533   // excluding current track
 534   // ---- beginning --- method from Andrea D 28.05.2010
 535   AliVertexerTracks *vertexer = new AliVertexerTracks(magneticField);
 536   vertexer->SetITSMode();
 537   vertexer->SetMinClusters(fNclustersITS);
 538   Int_t skipped[2];
 539   skipped[0] = (Int_t)track->GetID();
 540   vertexer->SetSkipTracks(1,skipped);
 541   AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent);
 542   delete vertexer; vertexer = NULL;
 543   if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return;
 544   // -- ending --- method from Andrea D 28.05.2010
 545
 546   Double_t dz[2];   // error of dca in cm
 547   Double_t covardz[3];
 548   if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
 549
 550   Double_t dcaXY = dz[0]*1.0e4;  // conv dca in cm to dca in micron
 551   Double_t dcaZ = dz[1]*1.0e4;
 552
 553   if(fDeDebugLevel>=10) printf("others----->electron %d has mother %d is from %.1f\n", eleLabel, ElectronFromSource(stack, eleLabel), var[0]);
 554
 555   if(TMath::Abs(dcaXY)<1000) var[1] = Int_t(dcaXY)/50;   // larger than 1mm should go to the last bin
 556   else
 557     var[1] = ((dcaXY>0)?1:-1)*20;
 558
 559   if(TMath::Abs(dcaZ)<1000) var[2] = Int_t(dcaZ)/100;
 560   else
 561     var[2] = ((dcaZ>0)?1:-1)*10;
 562
 563   // calculate dca using AliKFParticle class------------------------------------------------------------------
 564   Float_t  kfDcaXY = 0;
 565   Int_t trkID = track->GetID();
 566   AliKFParticle::SetField(magneticField);
 567   AliKFParticle kfParticle(*track, particle->GetPdgCode());
 568   // prepare kfprimary vertex
 569   AliKFVertex kfESDprimary;
 570   // Reconstruct Primary Vertex (with ESD tracks)
 571   Int_t n=primVtx->GetNIndices();
 572   if (n>0 && primVtx->GetStatus()){
 573     kfESDprimary = AliKFVertex(*primVtx);
 574     UShort_t *priIndex = primVtx->GetIndices();
 575     for (Int_t i=0;i<n;i++){
 576       Int_t idx = Int_t(priIndex[i]);
 577       if (idx == trkID){
 578         kfESDprimary -= kfParticle;
 579         kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4;
 580       }  // remove current track from this calculation
 581     }  // loop over all primary vertex contributors
 582   }
 583   // end of KF dca
 584
 585   if(TMath::Abs(kfDcaXY)<1000)
 586     var[3] = Int_t(kfDcaXY)/50;
 587   else
 588     var[3] = (kfDcaXY>0?1:-1)*20;
 589
 590   var[4] = pt;  // pt
 591   var[5] = eta; // eta
 592   var[6] = phi; // phi
 593   var[7] = mcR;
 594   var[8] = mcStatus;
 595
 596   (static_cast<THnSparseF *>(fDeOutputList->At(kEsdElectron)))->Fill(var);
 597
 598 }
 599
 600 //__________________________________________________________
 601 void AliHFEdisplacedElectrons::FillDataOutput(const AliESDEvent * const fESDEvent, AliESDtrack* const esdTrack)
 602 {
 603
 604   // this is pure data, without MC information at all
 605   // fill output: with HFE pid selection of electrons after all track quality cuts
 606
 607   if(!esdTrack) return;
 608   AliESDtrack *track = esdTrack;
 609
 610   Double_t pt   = track->Pt();
 611   Double_t eta = track->Eta();
 612   Double_t phi = track->Phi();
 613
 614   // obtain impact parameters in xy and y
 615   const AliESDVertex *primVtx = fESDEvent->GetPrimaryVertex();
 616
 617   Float_t magneticField = 5;  // initialized as 5kG
 618   magneticField = fESDEvent->GetMagneticField();  // in kG
 619   Double_t beampiperadius=3.;
 620
 621
 622   const Int_t nvarData = 6;
 623   Double_t varData[nvarData];
 624
 625   //
 626   // excluding current track
 627   //
 628
 629   //------ beginning --- method from Andrea D 28.05.2010
 630   AliVertexerTracks *vertexer = new AliVertexerTracks(fESDEvent->GetMagneticField());
 631   vertexer->SetITSMode();
 632   vertexer->SetMinClusters(fNclustersITS);
 633   Int_t skipped[2];
 634   skipped[0] = (Int_t)track->GetID();
 635   vertexer->SetSkipTracks(1,skipped);
 636   AliESDVertex *vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESDEvent);
 637   delete vertexer; vertexer = NULL;
 638   if(vtxESDSkip->GetNContributors()<fMinNprimVtxContributor) return;
 639   //------ ending --- method from Andrea D 28.05.2010
 640
 641   Double_t dz[2];   // error of dca in cm
 642   Double_t covardz[3];
 643
 644   if(!track->PropagateToDCA(vtxESDSkip,magneticField, beampiperadius, dz, covardz)) return; // protection
 645
 646   Double_t dcaXY = dz[0]*1.0e4;  // conv dca in cm to dca in micron
 647   Double_t dcaZ = dz[1]*1.0e4;
 648
 649   if(TMath::Abs(dcaXY)<1000) varData[0] = Int_t(dcaXY)/50;   // larger than 1mm should go to the last bin
 650   else
 651     varData[0] = (dcaXY>0?1:-1)*20;
 652   if(TMath::Abs(dcaZ)<1000) varData[1] = Int_t(dcaZ)/100;
 653   else
 654     varData[1] = (dcaZ>0?1:-1)*10;
 655
 656
 657   // calculate dca using AliKFParticle class------------------------------------------------------------------
 658   Float_t  kfDcaXY = 0;
 659   Int_t trkID = track->GetID();
 660   AliKFParticle::SetField(magneticField);
 661   AliKFParticle kfParticle(*track, -11*track->Charge());
 662   // prepare kfprimary vertex
 663   AliKFVertex kfESDprimary;
 664   // Reconstruct Primary Vertex (with ESD tracks)
 665   Int_t n=primVtx->GetNIndices();
 666   if (n>0 && primVtx->GetStatus()){
 667     kfESDprimary = AliKFVertex(*primVtx);
 668     UShort_t *priIndex = primVtx->GetIndices();
 669     for (Int_t i=0;i<n;i++){
 670       Int_t idx = Int_t(priIndex[i]);
 671       if (idx == trkID){
 672           kfESDprimary -= kfParticle;
 673           kfDcaXY = kfParticle.GetDistanceFromVertexXY(kfESDprimary)*1e4;
 674       }  // remove current track from this calculation
 675     }  // loop over all primary vertex contributors
 676   }
 677   if(TMath::Abs(kfDcaXY)<1000)
 678     varData[2] = Int_t(kfDcaXY)/50;
 679   else
 680     varData[2] = ((kfDcaXY)>0?1:-1)*20;
 681   // end of KF dca
 682
 683   varData[3] = pt; // pt
 684   varData[4] = eta; //eta
 685   varData[5] = phi; // phi
 686
 687   (static_cast<THnSparseF *>(fDeOutputList->At(kDataElectron)))->Fill(varData);
 688
 689 }
 690
 691 //__________________________________________________________
 692 Int_t AliHFEdisplacedElectrons::ElectronFromSource(AliStack * const stack, Int_t label) const
 693 {
 694
 695   //
 696   //  return electron source label via electron label
 697   //
 698
 699   // kEleSource is supposed to be either photon conv, direct photon conv, pi0, eta, beauty --> 0, 1, 2, 3, 4
 700
 701   Int_t eleLabel = label;
 702
 703   if(eleLabel<0 || eleLabel>stack->GetNtrack()) return -1;
 704
 705   TParticle *particle = stack->Particle(eleLabel);
 706   if(!particle) return -1;
 707
 708   Int_t motherLabel = particle->GetFirstMother();
 709   if(motherLabel<0 || motherLabel>stack->GetNtrack()) return -1;
 710   TParticle *motherPart = stack->Particle(motherLabel);
 711   if(!motherPart) return -1;
 712
 713   Int_t pdgCode = TMath::Abs(motherPart->GetPdgCode());
 714
 715   if(pdgCode==kPDGelectron) {
 716     if(fDeDebugLevel>=10) printf("particle label: %d...(motherLabel=%d : motherPdg=%d)  grandmother's pdg code was returned...%d \n",
 717                                label, motherLabel, pdgCode, ElectronFromSource(stack, motherLabel));
 718     return ElectronFromSource(stack, motherLabel);
 719   }
 720
 721   return pdgCode;
 722
 723 }
 724
 725 //__________________________________________________________
 726 Int_t AliHFEdisplacedElectrons::ElectronFromCharm(AliStack * const stack, Int_t label) const
 727 {
 728   //
 729   //  separate electron: kEleC from c->eX, kEleBC from b->c->eX
 730   //
 731
 732   Int_t motherLabel = label;
 733   TParticle *motherParticle = stack->Particle(motherLabel);  // mother part
 734   if(!motherParticle) return -1;
 735   Int_t gMotherLabel = motherParticle->GetFirstMother();  // grand mother
 736   if(gMotherLabel<0 || gMotherLabel>stack->GetNtrack()) return -1;
 737
 738   TParticle *gMotherPart = stack->Particle(gMotherLabel);
 739   if(!gMotherPart) return -1;
 740
 741   Int_t pdgCode = gMotherPart->GetPdgCode();
 742   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||    // for special intermediate meson states: like 100553
 743      TMath::Abs(pdgCode)/1000==kPDGbeauty || TMath::Abs(pdgCode)/100==kPDGbeauty || TMath::Abs(pdgCode)==kPDGbeauty) {
 744     if(fDeDebugLevel>=10)  printf("this electron label %d is from mother %d, and finally from %d\n", label, ElectronFromSource(stack, label), pdgCode);
 745     return kEleBC;
 746   }  // for sure it is from BC
 747
 748   else
 749     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||  // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 750        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 751        TMath::Abs(pdgCode)/100==kPDGcharm ||
 752        TMath::Abs(pdgCode)==kPDGcharm){
 753
 754       if(CharmFromBeauty(stack, gMotherLabel)!=-1) return kEleBC;
 755       else return kEleC;
 756
 757     }
 758
 759     else
 760       return -1;
 761
 762 }
 763
 764 //__________________________________________________________
 765 Int_t AliHFEdisplacedElectrons::CharmFromBeauty(AliStack * const stack, Int_t hfLabel) const
 766 {
 767
 768   //
 769   //  check if charm meson/hadron is from beauty decay
 770   //  -1, not from beauty
 771   //  return the label of the mother of this charm hadron
 772   //
 773
 774   Int_t charmLabel = hfLabel;
 775
 776   //  AliStack *stack = fMC->Stack();
 777   if(charmLabel<0 || charmLabel>stack->GetNtrack()) return -1;
 778   TParticle *particle = stack->Particle(charmLabel);
 779   if(!particle) return -1;
 780
 781   Int_t motherCharmLabel = particle->GetFirstMother();
 782   if(motherCharmLabel<0 || motherCharmLabel>stack->GetNtrack()) return -1;
 783
 784   TParticle *motherCharmPart = stack->Particle(motherCharmLabel);
 785   if(!motherCharmPart) return -1;
 786
 787   Int_t pdgCode = motherCharmPart->GetPdgCode();
 788
 789   if(TMath::Abs(pdgCode%10000)/100==kPDGbeauty ||   // for special intermediate meson states: like 100553
 790      TMath::Abs(pdgCode)/1000==kPDGbeauty ||
 791      TMath::Abs(pdgCode)/100==kPDGbeauty ||
 792      TMath::Abs(pdgCode)==kPDGbeauty)
 793     return motherCharmLabel;
 794   else
 795     if(TMath::Abs(pdgCode%10000)/100==kPDGcharm ||    // for special intermediate meson states: like 100443, 10443, 204*3 (*=1, 2, 3, 4)
 796        TMath::Abs(pdgCode)/1000==kPDGcharm ||
 797        TMath::Abs(pdgCode)/100==kPDGcharm)
 798       return CharmFromBeauty(stack, motherCharmLabel);   // loop over to see if charm is from beauty -- if yes, return the label of mother of the charm
 799
 800     else
 801       return -1;
 802 }
 803
 804 //__________________________________________________________
 805 Int_t AliHFEdisplacedElectrons::ElePhotonDirect(AliStack * const stack, Int_t label) const
 806 {
 807   //
 808   // electron is from photon, and check if this photon is direct
 809   //
 810
 811   Int_t eleLabel = label;
 812
 813   TParticle *particle = stack->Particle(eleLabel);
 814   if(particle->GetFirstMother()<0 || particle->GetFirstMother()>stack->GetNtrack())
 815     return -1;
 816
 817   TParticle *photonPart = stack->Particle(particle->GetFirstMother());
 818   if(!photonPart)
 819     return -1;
 820   Int_t motherPhotonLabel = photonPart->GetFirstMother();
 821   if(motherPhotonLabel<0 || motherPhotonLabel>stack->GetNtrack())
 822     return -1;
 823
 824   TParticle *motherPhotonPart = stack->Particle(motherPhotonLabel);
 825   if(!motherPhotonPart)
 826     return -1;
 827
 828   Int_t pdgMotherPhoton = motherPhotonPart->GetPdgCode();
 829   if(TMath::Abs(pdgMotherPhoton)<=10 || TMath::Abs(pdgMotherPhoton)==21)
 830     return 1;
 831
 832   else
 833     return -1;
 834
 835
 836 }
 837
 838 //__________________________________________________________
 839 Int_t AliHFEdisplacedElectrons::GetMCpid(AliStack* const stack, Int_t partLabel) const
 840 {
 841
 842   //
 843   // Simply pdg code
 844   //
 845
 846   Int_t label = partLabel;
 847 //   AliStack* stack = fMC->Stack();
 848   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 849
 850   // MC Information
 851   TParticle * particle = stack->Particle(label);
 852   if(!particle) return -1;
 853   Int_t pdg = particle->GetPdgCode();
 854
 855   return pdg;
 856
 857 }
 858
 859 //__________________________________________________________
 860 Int_t AliHFEdisplacedElectrons::GetMotherLabel(AliStack *const stack, Int_t eleLabel) const
 861 {
 862   //
 863   // Simply label of mother
 864   //
 865
 866
 867   Int_t label = eleLabel;
 868   //  AliStack* stack = fMC->Stack();
 869   if((label < 0) || (label >= stack->GetNtrack())) return -1;
 870
 871   // MC Information
 872   TParticle * particle = stack->Particle(label);
 873   if(!particle) return -1;
 874
 875   return particle->GetFirstMother();
 876
 877 }
 878
 879
 880 //__________________________________________________________
 881 Float_t AliHFEdisplacedElectrons::GetRapidity(TParticle *part) const
 882 {
 883   // return rapidity
 884
 885   Float_t rapidity = -999;
 886   if((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)
 887     rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 888
 889   return rapidity;
 890 }
 891
 892
 893 //__________________________________________________________
 894 Float_t AliHFEdisplacedElectrons::GetTrackRapidity(AliESDtrack * const track) const
 895 {
 896   // return rapidity of electron
 897
 898   Float_t px = track->Px();
 899   Float_t py = track->Py();
 900   Float_t pz = track->Pz();
 901
 902   // electron mass 0.00051099906 GeV/c2
 903   TParticlePDG* electron = TDatabasePDG::Instance()->GetParticle(kPDGelectron);
 904   Double_t mass = electron->Mass();
 905
 906   Float_t en = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);
 907
 908   Float_t rapidity = -999;
 909   if((en - pz)*(en + pz)>0)
 910     rapidity = 0.5*(TMath::Log((en - pz)*(en + pz)));
 911
 912   return rapidity;
 913 }
 914
 915 //__________________________________________________________
 916 Int_t AliHFEdisplacedElectrons::CheckCharm(AliStack * const stack, Int_t eleLabel)
 917 {
 918   // checks the genealogy of the ele from charm for a beauty
 919   // this method needs the stack and the label of the electron
 920   // warning: it assumes that the ele comes from a charm
 921
 922   TParticle *particle = stack->Particle(eleLabel);
 923   Int_t label = particle->GetFirstMother();
 924
 925   while(label>0 && label<(stack->GetNtrack())){
 926     particle = stack->Particle(label);
 927     if(IsB(TMath::Abs(particle->GetPdgCode()))) return kEleBC;
 928     label = particle->GetFirstMother();
 929   }
 930
 931   return kEleC;
 932 }
 933
 934 //__________________________________________________________
 935 Bool_t AliHFEdisplacedElectrons::IsB(Int_t pdg) const
 936 {
 937   // check if the pdg is that of a beauty particle
 938
 939   if((TMath::Abs(pdg)%10000)/100==kPDGbeauty ||
 940      TMath::Abs(pdg)/1000==kPDGbeauty ||
 941      TMath::Abs(pdg)/100==kPDGbeauty ||
 942      TMath::Abs(pdg)==kPDGbeauty) return kTRUE;
 943   else return kFALSE;
 944 }
 945
 946 //__________________________________________________________
 947 Bool_t AliHFEdisplacedElectrons::IsC(Int_t pdg) const
 948 {
 949   // check if the pdg is that of a charmed particle
 950
 951   if((TMath::Abs(pdg)%10000)/100==kPDGcharm ||
 952      TMath::Abs(pdg)/1000==kPDGcharm ||
 953      TMath::Abs(pdg)/100==kPDGcharm ||
 954      TMath::Abs(pdg)==kPDGcharm) return kTRUE;
 955   else return kFALSE;
 956 }