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