PWG4/PartCorrDep/AliAnaElectron.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 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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for the electron identification.
  20 // Clusters from EMCAL matched to tracks
  21 // and kept in the AOD. Few histograms produced.
  22 //
  23 // -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26 // --- ROOT system ---
  27 #include <TH2F.h>
  28 #include <TParticle.h>
  29 #include <TNtuple.h>
  30 #include <TClonesArray.h>
  31 #include <TObjString.h>
  32 //#include <Riostream.h>
  33
  34 // --- Analysis system ---
  35 #include "AliAnaElectron.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliMCAnalysisUtils.h"
  38 #include "AliAODCaloCluster.h"
  39 #include "AliFidutialCut.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODPid.h"
  42 #include "AliCaloPID.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliStack.h"
  45 #include "AliExternalTrackParam.h"
  46 #include "AliESDv0.h"
  47
  48 ClassImp(AliAnaElectron)
  49
  50 //____________________________________________________________________________
  51 AliAnaElectron::AliAnaElectron()
  52 : AliAnaPartCorrBaseClass(),fCalorimeter(""),
  53   fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),
  54   fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),
  55   fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),
  56   fWriteNtuple(kFALSE),
  57   //matching checks
  58   fEleNtuple(0),
  59   fh1pOverE(0),fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),
  60   fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),
  61   fh2TrackPVsClusterE(0),fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),
  62   //Photonic electron checks
  63   fh1OpeningAngle(0),fh1MinvPhoton(0),
  64   //reco
  65   fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
  66   fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),
  67   fhPtPE(0),fhPhiPE(0),fhEtaPE(0),
  68   fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
  69   fhPtBottom(0),fhPhiBottom(0),fhEtaBottom(0),
  70   fhPtCharm(0),fhPhiCharm(0),fhEtaCharm(0),
  71   fhPtCFromB(0),fhPhiCFromB(0),fhEtaCFromB(0),
  72   fhPtDalitz(0),fhPhiDalitz(0),fhEtaDalitz(0),
  73   fhPtWDecay(0),fhPhiWDecay(0),fhEtaWDecay(0),
  74   fhPtZDecay(0),fhPhiZDecay(0),fhEtaZDecay(0),
  75   fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
  76   fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
  77   fhPtHadron(0),fhPtEleTrkDet(0),
  78   //B-tagging
  79   fhBtagCut1(0),fhBtagCut2(0),fhBtagCut3(0),fhBtagQA1(0),fhBtagQA2(0),
  80   //MC
  81   fMCEleNtuple(0),fhPtMCHadron(0)
  82 {
  83   //default ctor
  84
  85   //Initialize parameters
  86   InitParameters();
  87
  88 }
  89
  90 //____________________________________________________________________________
  91 AliAnaElectron::AliAnaElectron(const AliAnaElectron & g)
  92   : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),
  93    fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),fResidualCut(g.fResidualCut),
  94    fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),
  95   fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),
  96    fWriteNtuple(g.fWriteNtuple),
  97    //matching checks
  98    fEleNtuple(g.fEleNtuple),
  99    fh1pOverE(g.fh1pOverE),fh1dR(g.fh1dR),
 100    fh2EledEdx(g.fh2EledEdx),fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),
 101    fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),
 102    fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),
 103    fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),
 104    //Photonic electron checks
 105    fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),
 106    //reco
 107    fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron),
 108    fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE),
 109    fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE),
 110    fhPtConversion(g.fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
 111    fhPtBottom(g.fhPtBottom),fhPhiBottom(g.fhPhiBottom),fhEtaBottom(g.fhEtaBottom),
 112    fhPtCharm(g.fhPtCharm),fhPhiCharm(g.fhPhiCharm),fhEtaCharm(g.fhEtaCharm),
 113    fhPtCFromB(g.fhPtCFromB),fhPhiCFromB(g.fhPhiCFromB),fhEtaCFromB(g.fhEtaCFromB),
 114    fhPtDalitz(g.fhPtDalitz),fhPhiDalitz(g.fhPhiDalitz),fhEtaDalitz(g.fhEtaDalitz),
 115    fhPtWDecay(g.fhPtWDecay),fhPhiWDecay(g.fhPhiWDecay),fhEtaWDecay(g.fhEtaWDecay),
 116    fhPtZDecay(g.fhPtZDecay),fhPhiZDecay(g.fhPhiZDecay),fhEtaZDecay(g.fhEtaZDecay),
 117    fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
 118    fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
 119    fhPtHadron(g.fhPtHadron),fhPtEleTrkDet(g.fhPtEleTrkDet),
 120    //B-tagging
 121    fhBtagCut1(g.fhBtagCut1),fhBtagCut2(g.fhBtagCut2),fhBtagCut3(g.fhBtagCut3),
 122    fhBtagQA1(g.fhBtagQA1),fhBtagQA2(g.fhBtagQA2),
 123    //MC
 124     fMCEleNtuple(g.fMCEleNtuple),fhPtMCHadron(g.fhPtMCHadron)
 125 {
 126   // cpy ctor
 127
 128 }
 129
 130 //_________________________________________________________________________
 131 AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)
 132 {
 133   // assignment operator
 134
 135   if(&g == this) return *this;
 136   fCalorimeter = g.fCalorimeter;
 137   fpOverEmin = g.fpOverEmin;
 138   fpOverEmax = g.fpOverEmax;
 139   fResidualCut = g.fResidualCut;
 140   fDrCut = g.fDrCut;
 141   fPairDcaCut = g.fPairDcaCut;
 142   fDecayLenCut = g.fDecayLenCut;
 143   fImpactCut = g.fImpactCut;
 144   fAssocPtCut = g.fAssocPtCut;
 145   fMassCut = g.fMassCut;
 146   fSdcaCut = g.fSdcaCut;
 147   fITSCut = g.fITSCut;
 148   fWriteNtuple = g.fWriteNtuple;
 149   fEleNtuple = g.fEleNtuple;
 150   fh1pOverE = g.fh1pOverE;
 151   fh1dR = g.fh1dR;
 152   fh2EledEdx = g.fh2EledEdx;
 153   fh2MatchdEdx = g.fh2MatchdEdx;
 154   fh2dEtadPhi = g.fh2dEtadPhi;
 155   fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;
 156   fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;
 157   fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;
 158   fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;
 159   fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;
 160   fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;
 161   fh1OpeningAngle = g.fh1OpeningAngle;
 162   fh1MinvPhoton = g.fh1MinvPhoton;
 163   fhPtElectron = g.fhPtElectron;
 164   fhPhiElectron = g.fhPhiElectron;
 165   fhEtaElectron = g.fhEtaElectron;
 166   fhPtNPE = g.fhPtNPE;
 167   fhPhiNPE = g.fhPhiNPE;
 168   fhEtaNPE = g.fhEtaNPE;
 169   fhPtPE = g.fhPtPE;
 170   fhPhiPE = g.fhPhiPE;
 171   fhEtaPE = g.fhEtaPE;
 172   fhPtConversion = g.fhPtConversion;
 173   fhPhiConversion = g.fhPhiConversion;
 174   fhEtaConversion = g.fhEtaConversion;
 175   fhPtBottom = g.fhPtBottom;
 176   fhPhiBottom = g.fhPhiBottom;
 177   fhEtaBottom = g.fhEtaBottom;
 178   fhPtCharm = g.fhPtCharm;
 179   fhPhiCharm = g.fhPhiCharm;
 180   fhEtaCharm = g.fhEtaCharm;
 181   fhPtCFromB = g.fhPtCFromB;
 182   fhPhiCFromB = g.fhPhiCFromB;
 183   fhEtaCFromB = g.fhEtaCFromB;
 184   fhPtDalitz = g.fhPtDalitz;
 185   fhPhiDalitz = g.fhPhiDalitz;
 186   fhEtaDalitz = g.fhEtaDalitz;
 187   fhPtWDecay = g.fhPtWDecay;
 188   fhPhiWDecay = g.fhPhiWDecay;
 189   fhEtaWDecay = g.fhEtaWDecay;
 190   fhPtZDecay = g.fhPtZDecay;
 191   fhPhiZDecay = g.fhPhiZDecay;
 192   fhEtaZDecay = g.fhEtaZDecay;
 193   fhPtPrompt = g.fhPtPrompt;
 194   fhPhiPrompt = g.fhPhiPrompt;
 195   fhEtaPrompt = g.fhEtaPrompt;
 196   fhPtUnknown = g.fhPtUnknown;
 197   fhPhiUnknown = g.fhPhiUnknown;
 198   fhEtaUnknown = g.fhEtaUnknown;
 199
 200   fhPtHadron = g.fhPtHadron;
 201   fhPtEleTrkDet = g.fhPtEleTrkDet;
 202
 203   //B-tagging
 204   fhBtagCut1 = g.fhBtagCut1;
 205   fhBtagCut2 = g.fhBtagCut2;
 206   fhBtagCut3 = g.fhBtagCut3;
 207   fhBtagQA1 = g.fhBtagQA1;
 208   fhBtagQA2 = g.fhBtagQA2;
 209
 210   fMCEleNtuple = g.fMCEleNtuple;
 211   fhPtMCHadron = g.fhPtMCHadron;
 212
 213   return *this;
 214
 215 }
 216
 217 //____________________________________________________________________________
 218 AliAnaElectron::~AliAnaElectron()
 219 {
 220   //dtor
 221
 222 }
 223
 224
 225 //________________________________________________________________________
 226 TList *  AliAnaElectron::GetCreateOutputObjects()
 227 {
 228   // Create histograms to be saved in output file and
 229   // store them in outputContainer
 230   TList * outputContainer = new TList() ;
 231   outputContainer->SetName("ElectronHistos") ;
 232
 233   //created ele ntuple for further analysis
 234   if(fWriteNtuple) {
 235       fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");
 236     outputContainer->Add(fEleNtuple) ;
 237   }
 238
 239   Int_t nptbins  = GetHistoNPtBins();
 240   Int_t nphibins = GetHistoNPhiBins();
 241   Int_t netabins = GetHistoNEtaBins();
 242   Float_t ptmax  = GetHistoPtMax();
 243   Float_t phimax = GetHistoPhiMax();
 244   Float_t etamax = GetHistoEtaMax();
 245   Float_t ptmin  = GetHistoPtMin();
 246   Float_t phimin = GetHistoPhiMin();
 247   Float_t etamin = GetHistoEtaMin();
 248
 249   fh1pOverE = new TH1F("h1pOverE","EMCAL-TRACK matches p/E",100,0.,10.);
 250   fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi());
 251   fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.);
 252   fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.);
 253   fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
 254   fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
 255   fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
 256
 257   fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 258   fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 259   fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax);
 260   fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax);
 261
 262   outputContainer->Add(fh1pOverE) ;
 263   outputContainer->Add(fh1dR) ;
 264   outputContainer->Add(fh2EledEdx) ;
 265   outputContainer->Add(fh2MatchdEdx) ;
 266   outputContainer->Add(fh2dEtadPhi) ;
 267   outputContainer->Add(fh2dEtadPhiMatched) ;
 268   outputContainer->Add(fh2dEtadPhiUnmatched) ;
 269   outputContainer->Add(fh2TrackPVsClusterE) ;
 270   outputContainer->Add(fh2TrackPtVsClusterE) ;
 271   outputContainer->Add(fh2TrackPhiVsClusterPhi) ;
 272   outputContainer->Add(fh2TrackEtaVsClusterEta) ;
 273
 274   //photonic electron checks
 275   fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between e+e- pairs",100,0.,TMath::Pi());
 276   fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of e+e- pairs",200,0.,2.);
 277
 278   outputContainer->Add(fh1OpeningAngle);
 279   outputContainer->Add(fh1MinvPhoton);
 280
 281   fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);
 282   fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 283   fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 284   fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);
 285   fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 286   fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 287   fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);
 288   fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 289   fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 290
 291   outputContainer->Add(fhPtElectron) ;
 292   outputContainer->Add(fhPhiElectron) ;
 293   outputContainer->Add(fhEtaElectron) ;
 294   outputContainer->Add(fhPtNPE) ;
 295   outputContainer->Add(fhPhiNPE) ;
 296   outputContainer->Add(fhEtaNPE) ;
 297   outputContainer->Add(fhPtPE) ;
 298   outputContainer->Add(fhPhiPE) ;
 299   outputContainer->Add(fhEtaPE) ;
 300
 301   fhPtHadron = new TH1F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);
 302   fhPtEleTrkDet = new TH1F("hPtEleTrkDet","Electrons identified by tracking detectors w/in EMCAL acceptance",nptbins,ptmin,ptmax);
 303
 304   outputContainer->Add(fhPtHadron);
 305   outputContainer->Add(fhPtEleTrkDet);
 306
 307   //B-tagging
 308   fhBtagCut1 = new TH2F("hbtag_cut1","B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);
 309   fhBtagCut2 = new TH2F("hbtag_cut2","B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);
 310   fhBtagCut3 = new TH2F("hbtag_cut3","B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);
 311   fhBtagQA1  = new TH2F("hbtag_qa1" ,"B-tag QA: pairDCA vs length", 100,0,0.2 ,100,0,1.0);
 312   fhBtagQA2  = new TH2F("hbtag_qa2" ,"B-tag QA: signDCA vs mass"  , 200,-0.5,0.5 ,100,0,10);
 313
 314   outputContainer->Add(fhBtagCut1) ;
 315   outputContainer->Add(fhBtagCut2) ;
 316   outputContainer->Add(fhBtagCut3) ;
 317   outputContainer->Add(fhBtagQA1) ;
 318   outputContainer->Add(fhBtagQA2) ;
 319
 320   if(IsDataMC()){
 321
 322     fhPtConversion = new TH1F("hPtConversion","Conversion electron pT",nptbins,ptmin,ptmax);
 323     fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 324     fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 325     fhPtBottom = new TH1F("hPtBottom","Bottom electron pT",nptbins,ptmin,ptmax);
 326     fhPhiBottom = new TH2F("hPhiBottom","Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 327     fhEtaBottom = new TH2F("hEtaBottom","Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 328     fhPtCharm = new TH1F("hPtCharm","Charm electron pT",nptbins,ptmin,ptmax);
 329     fhPhiCharm = new TH2F("hPhiCharm","Charm Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 330     fhEtaCharm = new TH2F("hEtaCharm","Charm Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 331     fhPtCFromB = new TH1F("hPtCFromB","Charm from Bottom electron pT",nptbins,ptmin,ptmax);
 332     fhPhiCFromB = new TH2F("hPhiCFromB","Charm from Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 333     fhEtaCFromB = new TH2F("hEtaCFromB","Charm from Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 334     fhPtDalitz = new TH1F("hPtDalitz","Dalitz electron pT",nptbins,ptmin,ptmax);
 335     fhPhiDalitz = new TH2F("hPhiDalitz","Dalitz Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 336     fhEtaDalitz = new TH2F("hEtaDalitz","Dalitz Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 337     fhPtWDecay = new TH1F("hPtWDecay","W-boson Electron pT",nptbins,ptmin,ptmax);
 338     fhPhiWDecay = new TH2F("hPhiWDecay","W-boson electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 339     fhEtaWDecay = new TH2F("hEtaWDecay","W-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 340     fhPtZDecay = new TH1F("hPtZDecay","Z-boson electron pT",nptbins,ptmin,ptmax);
 341     fhPhiZDecay = new TH2F("hPhiZDecay","Z-boson Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 342     fhEtaZDecay = new TH2F("hEtaZDecay","Z-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 343     fhPtPrompt = new TH1F("hPtPrompt","Prompt electron pT",nptbins,ptmin,ptmax);
 344     fhPhiPrompt = new TH2F("hPhiPrompt","Prompt Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 345     fhEtaPrompt = new TH2F("hEtaPrompt","Prompt Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 346     fhPtUnknown = new TH1F("hPtUnknown","Unknown electron pT",nptbins,ptmin,ptmax);
 347     fhPhiUnknown = new TH2F("hPhiUnknown","Unknown Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 348     fhEtaUnknown = new TH2F("hEtaUnknown","Unknown Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 349
 350     outputContainer->Add(fhPtConversion);
 351     outputContainer->Add(fhPhiConversion);
 352     outputContainer->Add(fhEtaConversion);
 353     outputContainer->Add(fhPtBottom);
 354     outputContainer->Add(fhPhiBottom);
 355     outputContainer->Add(fhEtaBottom);
 356     outputContainer->Add(fhPtCharm);
 357     outputContainer->Add(fhPhiCharm);
 358     outputContainer->Add(fhEtaCharm);
 359     outputContainer->Add(fhPtCFromB);
 360     outputContainer->Add(fhPhiCFromB);
 361     outputContainer->Add(fhEtaCFromB);
 362     outputContainer->Add(fhPtDalitz);
 363     outputContainer->Add(fhPhiDalitz);
 364     outputContainer->Add(fhEtaDalitz);
 365     outputContainer->Add(fhPtWDecay);
 366     outputContainer->Add(fhPhiWDecay);
 367     outputContainer->Add(fhEtaWDecay);
 368     outputContainer->Add(fhPtZDecay);
 369     outputContainer->Add(fhPhiZDecay);
 370     outputContainer->Add(fhEtaZDecay);
 371     outputContainer->Add(fhPtPrompt);
 372     outputContainer->Add(fhPhiPrompt);
 373     outputContainer->Add(fhEtaPrompt);
 374     outputContainer->Add(fhPtUnknown);
 375     outputContainer->Add(fhPhiUnknown);
 376     outputContainer->Add(fhEtaUnknown);
 377
 378     //created ele ntuple for further analysis
 379     if(fWriteNtuple) {
 380       fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");
 381       outputContainer->Add(fMCEleNtuple) ;
 382     }
 383
 384     fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);
 385
 386     outputContainer->Add(fhPtMCHadron);
 387
 388   }//Histos with MC
 389
 390   //Save parameters used for analysis
 391   TString parList ; //this will be list of parameters used for this analysis.
 392   char onePar[255] ;
 393
 394   sprintf(onePar,"--- AliAnaElectron ---\n") ;
 395   parList+=onePar ;
 396   sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ;
 397   parList+=onePar ;
 398   sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ;
 399   parList+=onePar ;
 400   sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ;
 401   parList+=onePar ;
 402   sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ;
 403   parList+=onePar ;
 404   sprintf(onePar,"---Btagging\n");
 405   parList+=onePar ;
 406   sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut);
 407   parList+=onePar ;
 408   sprintf(onePar,"min ITS-hits: %d\n",fITSCut);
 409   parList+=onePar ;
 410   sprintf(onePar,"max dR (e,h): %f\n",fDrCut);
 411   parList+=onePar ;
 412   sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut);
 413   parList+=onePar ;
 414   sprintf(onePar,"max decaylength: %f\n",fDecayLenCut);
 415   parList+=onePar ;
 416   sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut);
 417   parList+=onePar ;
 418
 419   //Get parameters set in base class.
 420   parList += GetBaseParametersList() ;
 421
 422   //Get parameters set in FidutialCut class (not available yet)
 423   //parlist += GetFidCut()->GetFidCutParametersList()
 424
 425   TObjString *oString= new TObjString(parList) ;
 426   outputContainer->Add(oString);
 427
 428   return outputContainer ;
 429
 430 }
 431
 432 //____________________________________________________________________________
 433 void AliAnaElectron::Init()
 434 {
 435
 436   //do some initialization
 437   if(fCalorimeter == "PHOS") {
 438     printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");
 439     fCalorimeter = "EMCAL";
 440   }
 441   if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
 442     printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");
 443     abort();
 444   }
 445
 446 }
 447
 448
 449 //____________________________________________________________________________
 450 void AliAnaElectron::InitParameters()
 451 {
 452
 453   //Initialize the parameters of the analysis.
 454   SetOutputAODClassName("AliAODPWG4Particle");
 455   SetOutputAODName("PWG4Particle");
 456
 457   AddToHistogramsName("AnaElectron_");
 458
 459   fCalorimeter = "EMCAL" ;
 460   fpOverEmin = 0.5;
 461   fpOverEmax = 1.5;
 462   fResidualCut = 0.02;
 463   //B-tagging
 464   fDrCut       = 1.0;
 465   fPairDcaCut  = 0.02;
 466   fDecayLenCut = 1.0;
 467   fImpactCut   = 0.5;
 468   fAssocPtCut  = 1.0;
 469   fMassCut     = 1.5;
 470   fSdcaCut     = 0.1;
 471   fITSCut      = 4;
 472
 473 }
 474
 475 //__________________________________________________________________
 476 void  AliAnaElectron::MakeAnalysisFillAOD()
 477 {
 478   //
 479   // Do analysis and fill aods with electron candidates
 480   // These AODs will be used to do subsequent histogram filling
 481   //
 482   // Also fill some QA histograms
 483   //
 484
 485   TObjArray *cl = new TObjArray();
 486
 487   Double_t bfield = 0.;
 488   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
 489
 490   //Select the calorimeter of the electron
 491   if(fCalorimeter != "EMCAL") {
 492     printf("This class not yet implemented for PHOS\n");
 493     abort();
 494   }
 495   cl = GetAODEMCAL();
 496
 497   ////////////////////////////////////////////////
 498   //Start from tracks and get associated clusters
 499   ////////////////////////////////////////////////
 500   if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
 501   Int_t ntracks = GetAODCTS()->GetEntriesFast();
 502   if(GetDebug() > 0)
 503     printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
 504
 505   //Unfortunately, AliAODTracks don't have associated EMCAL clusters.
 506   //we have to redo track-matching, I guess
 507   Int_t iCluster = -999;
 508   Int_t bt = 0; //counter for event b-tags
 509
 510   for (Int_t itrk =  0; itrk <  ntracks; itrk++) {////////////// track loop
 511     iCluster = -999; //start with no match
 512     AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;
 513     AliAODPid* pid = (AliAODPid*) track->GetDetPid();
 514
 515     Double_t emcpos[3];
 516     pid->GetEMCALPosition(emcpos);
 517     Double_t emcmom[3];
 518     pid->GetEMCALMomentum(emcmom);
 519
 520     TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);
 521     TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);
 522     Double_t tphi = pos.Phi();
 523     Double_t teta = pos.Eta();
 524     Double_t tmom = mom.Mag();
 525
 526     TLorentzVector mom2(mom,0.);
 527     Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter) ;
 528     if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, eta %2.2f in fidutial cut %d\n",track->Pt(), track->Phi(), track->Eta(), in);
 529     if(mom.Pt() > GetMinPt() && in) {
 530
 531       Double_t dEdx = pid->GetTPCsignal();
 532
 533       //NOTE:  As of 02-Sep-2009, the XYZAtDCA methods of AOD do not
 534       //work, but it is possible to get the position of a track at
 535       //closest approach to the vertex from the GetPosition method
 536       Double_t xyz[3];
 537       //track->XYZAtDCA(xyz);
 538       Bool_t isNotDCA = track->GetPosition(xyz);
 539       if(isNotDCA) printf("##Problem getting impact parameter!\n");
 540       //printf("\tTRACK POSITION AT DCA: %2.2f,%2.2f,%2.2f\n",xyz[0],xyz[1],xyz[2]);
 541       Double_t xy = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 542       Double_t z = xyz[2];
 543
 544       Int_t ntot = cl->GetEntriesFast();
 545       Double_t res = 999.;
 546       Double_t pOverE = -999.;
 547
 548       Int_t pidProb = track->GetMostProbablePID();
 549       if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) fhPtHadron->Fill(track->Pt());
 550       if(pidProb == AliAODTrack::kElectron) fhPtEleTrkDet->Fill(track->Pt());
 551
 552       Bool_t isElectron = kFALSE;
 553       //For tracks in EMCAL acceptance, pair them with all clusters
 554       //and fill the dEta vs dPhi for these pairs:
 555       for(Int_t iclus = 0; iclus < ntot; iclus++) {
 556         AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));
 557         if(!clus) continue;
 558
 559         Double_t x[3];
 560         clus->GetPosition(x);
 561         TVector3 cluspos(x[0],x[1],x[2]);
 562         Double_t deta = teta - cluspos.Eta();
 563         Double_t dphi = tphi - cluspos.Phi();
 564         if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
 565         if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
 566         fh2dEtadPhi->Fill(deta,dphi);
 567         fh2TrackPVsClusterE->Fill(clus->E(),track->P());
 568         fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());
 569         fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());
 570         fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());
 571
 572         res = sqrt(dphi*dphi + deta*deta);
 573         fh1dR->Fill(res);
 574
 575         /////////////////////////////////
 576         //Perform electron cut analysis//
 577         /////////////////////////////////
 578         //Good match
 579         if(res < fResidualCut) {
 580           fh2dEtadPhiMatched->Fill(deta,dphi);
 581           iCluster = iclus;
 582
 583           Int_t tmctag = -1;
 584           Int_t cmctag = -1;
 585
 586           if(IsDataMC()) {
 587             //Input from second AOD?
 588             Int_t input = 0;
 589             if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;
 590             tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);
 591
 592             //Do you want the cluster or the track label?
 593             input = 0;
 594             if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;
 595             cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(0),GetReader(),input);
 596           }
 597
 598           if(fWriteNtuple) {
 599             fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,pidProb,xy,z);
 600           }
 601
 602           fh2MatchdEdx->Fill(track->P(),dEdx);
 603
 604           Double_t energy = clus->E();
 605           if(energy > 0) pOverE = tmom/energy;
 606           fh1pOverE->Fill(pOverE);
 607
 608           Int_t mult = clus->GetNCells();
 609           if(mult < 2 &&  GetDebug() > 0) printf("Single digit cluster.\n");
 610
 611           //////////////////////////////
 612           //Electron cuts happen here!//
 613           //////////////////////////////
 614           if(pOverE > fpOverEmin && pOverE < fpOverEmax) isElectron = kTRUE;
 615         } else {
 616           fh2dEtadPhiUnmatched->Fill(deta,dphi);
 617         }
 618
 619       } //calocluster loop
 620
 621       ///////////////////////////
 622       //Fill AOD with electrons//
 623       ///////////////////////////
 624       if(isElectron) {
 625
 626         //B-tagging
 627         if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");
 628         Int_t btag = GetBtag(track); bt += btag;
 629
 630         fh2EledEdx->Fill(track->P(),dEdx);
 631
 632         Double_t eMass = 0.511/1000; //mass in GeV
 633         Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);
 634         AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);
 635         tr.SetLabel(track->GetLabel());
 636         tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters
 637         tr.SetTrackLabel(itrk,-1); //sets the indices of the original tracks
 638         tr.SetDetector(fCalorimeter);
 639         if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);
 640         //Make this preserve sign of particle
 641         if(track->Charge() < 0) tr.SetPdg(11); //electron is 11
 642         else  tr.SetPdg(-11); //positron is -11
 643         tr.SetBtag(btag);
 644
 645         //Play with the MC stack if available
 646         //Check origin of the candidates
 647         if(IsDataMC()){
 648
 649           //FIXME:  Need to re-think this for track-oriented analysis
 650           //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?
 651           tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));
 652
 653           if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());
 654         }//Work with stack also
 655
 656         AddAODParticle(tr);
 657
 658         if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());
 659       }//electron
 660     }//pt, fiducial selection
 661   }//track loop
 662
 663   //FIXME:  Should we also check from the calocluster side, just in
 664   //case?
 665
 666   if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");
 667   if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD()  End fill AODs \n");
 668
 669 }
 670
 671 //__________________________________________________________________
 672 void  AliAnaElectron::MakeAnalysisFillHistograms()
 673 {
 674   //Do analysis and fill histograms
 675
 676   AliStack * stack = 0x0;
 677   TParticle * primary = 0x0;
 678   TClonesArray * mcparticles0 = 0x0;
 679   TClonesArray * mcparticles1 = 0x0;
 680   AliAODMCParticle * aodprimary = 0x0;
 681
 682   if(IsDataMC()) {
 683     if(GetReader()->ReadStack()){
 684       stack =  GetMCStack() ;
 685
 686       if(!stack)
 687         printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");
 688
 689     }
 690     else if(GetReader()->ReadAODMCParticles()){
 691       //Get the list of MC particles
 692       mcparticles0 = GetReader()->GetAODMCParticles(0);
 693       if(!mcparticles0 && GetDebug() > 0)     {
 694         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 695       }
 696       if(GetReader()->GetSecondInputAODTree()){
 697         mcparticles1 = GetReader()->GetAODMCParticles(1);
 698         if(!mcparticles1 && GetDebug() > 0)     {
 699           printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
 700         }
 701       }
 702
 703     }
 704   }// is data and MC
 705
 706   //Loop on stored AOD electrons
 707   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 708   if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 709
 710   for(Int_t iaod = 0; iaod < naod ; iaod++){
 711     AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 712     Int_t pdg = ele->GetPdg();
 713
 714     if(GetDebug() > 3)
 715       printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;
 716
 717     if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue;
 718     if(ele->GetDetector() != fCalorimeter) continue;
 719
 720     if(GetDebug() > 1)
 721       printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;
 722
 723
 724     //Filter for photonic electrons based on opening angle and Minv
 725     //cuts, also fill histograms
 726     Bool_t photonic = kFALSE;
 727     photonic = IsItPhotonic(ele);
 728
 729     //Fill electron histograms
 730     Float_t ptele = ele->Pt();
 731     Float_t phiele = ele->Phi();
 732     Float_t etaele = ele->Eta();
 733
 734     fhPtElectron  ->Fill(ptele);
 735     fhPhiElectron ->Fill(ptele,phiele);
 736     fhEtaElectron ->Fill(ptele,etaele);
 737
 738     if(photonic) {
 739       fhPtPE->Fill(ptele);
 740       fhPhiPE->Fill(ptele,phiele);
 741       fhEtaPE->Fill(ptele,etaele);
 742     } else {
 743       fhPtNPE->Fill(ptele);
 744       fhPhiNPE->Fill(ptele,phiele);
 745       fhEtaNPE->Fill(ptele,etaele);
 746     }
 747
 748     if(IsDataMC()){
 749       Int_t tag = ele->GetTag();
 750       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)){
 751         fhPtConversion  ->Fill(ptele);
 752         fhPhiConversion ->Fill(ptele,phiele);
 753         fhEtaConversion ->Fill(ptele,etaele);
 754       }
 755       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB))
 756         {
 757           fhPtBottom  ->Fill(ptele);
 758           fhPhiBottom ->Fill(ptele,phiele);
 759           fhEtaBottom ->Fill(ptele,etaele);
 760         }
 761       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC))
 762         {
 763           fhPtCharm  ->Fill(ptele);
 764           fhPhiCharm ->Fill(ptele,phiele);
 765           fhEtaCharm ->Fill(ptele,etaele);
 766         }
 767       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB))
 768         {
 769           fhPtCFromB  ->Fill(ptele);
 770           fhPhiCFromB ->Fill(ptele,phiele);
 771           fhEtaCFromB ->Fill(ptele,etaele);
 772         }
 773       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 774         {
 775           fhPtDalitz  ->Fill(ptele);
 776           fhPhiDalitz ->Fill(ptele,phiele);
 777           fhEtaDalitz ->Fill(ptele,etaele);
 778         }
 779       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay))
 780         {
 781           fhPtWDecay  ->Fill(ptele);
 782           fhPhiWDecay ->Fill(ptele,phiele);
 783           fhEtaWDecay ->Fill(ptele,etaele);
 784         }
 785       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay))
 786         {
 787           fhPtZDecay  ->Fill(ptele);
 788           fhPhiZDecay ->Fill(ptele,phiele);
 789           fhEtaZDecay ->Fill(ptele,etaele);
 790         }
 791       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
 792         {
 793           fhPtPrompt  ->Fill(ptele);
 794           fhPhiPrompt ->Fill(ptele,phiele);
 795           fhEtaPrompt ->Fill(ptele,etaele);
 796         }
 797       else{
 798         fhPtUnknown  ->Fill(ptele);
 799         fhPhiUnknown ->Fill(ptele,phiele);
 800         fhEtaUnknown ->Fill(ptele,etaele);
 801       }
 802     }//Histograms with MC
 803
 804   }// aod loop
 805
 806   ////////////////////////////////////////////////////////
 807   //Fill histograms of pure MC kinematics from the stack//
 808   ////////////////////////////////////////////////////////
 809   if(IsDataMC()) {
 810     if(GetReader()->ReadStack()) {
 811       for(Int_t ipart = 0; ipart < stack->GetNtrack(); ipart++) {
 812         primary = stack->Particle(ipart);
 813         TLorentzVector mom;
 814         primary->Momentum(mom);
 815         Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter);
 816         if(primary->Pt() < GetMinPt()) continue;
 817         if(!in) continue;
 818
 819         Int_t pdgcode = primary->GetPdgCode();
 820         if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212)
 821           fhPtMCHadron->Fill(primary->Pt());
 822
 823         //we only care about electrons
 824         if(TMath::Abs(pdgcode) != 11) continue;
 825         //we only want TRACKABLE electrons (TPC 85-250cm)
 826         if(primary->R() > 200.) continue;
 827         //Ignore low pt electrons
 828         if(primary->Pt() < 0.2) continue;
 829
 830         //find out what the ancestry of this electron is
 831         Int_t mctag = -1;
 832         Int_t input = 0;
 833         mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);
 834
 835         //fill ntuple
 836         if(fWriteNtuple) {
 837           fMCEleNtuple->Fill(mctag,primary->Pt(),primary->Phi(),primary->Eta(),primary->Vx(),primary->Vy(),primary->Vz());
 838         }
 839
 840       }
 841
 842     } else if(GetReader()->ReadAODMCParticles()) {
 843       Int_t npart0 = mcparticles0->GetEntriesFast();
 844       Int_t npart1 = 0;
 845       if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();
 846       Int_t npart = npart0+npart1;
 847       for(Int_t ipart = 0; ipart < npart; ipart++) {
 848         if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);
 849         else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);
 850         if(!aodprimary) {
 851           printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", ipart);
 852           continue;
 853         }
 854
 855         Double_t mom[3] = {0.,0.,0.};
 856         aodprimary->PxPyPz(mom);
 857         TLorentzVector mom2(mom,0.);
 858         Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter);
 859         if(aodprimary->Pt() < GetMinPt()) continue;
 860         if(!in) continue;
 861
 862         Int_t pdgcode = aodprimary->GetPdgCode();
 863         if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212)
 864           fhPtMCHadron->Fill(aodprimary->Pt());
 865
 866         //we only care about electrons
 867         if(TMath::Abs(pdgcode) != 11) continue;
 868         //we only want TRACKABLE electrons (TPC 85-250cm)
 869         Double_t radius = TMath::Sqrt(aodprimary->Xv()*aodprimary->Xv() + aodprimary->Yv()*aodprimary->Yv());
 870         if(radius > 200.) continue;
 871
 872         if(aodprimary->Pt() < 0.2) continue;
 873
 874         /*
 875         if(aodprimary->GetDaughter(0) > 0) {
 876           Int_t dindex = aodprimary->GetDaughter(0);
 877           //printf("####AODMCparticle daughter index %d and flag value %d\n",dindex,aodprimary->GetFlag());
 878           Double_t dxv = 0.;
 879           Double_t dyv = 0.;
 880           if(ipart < npart0) {
 881             dxv =  ((AliAODMCParticle*)mcparticles0->At(dindex))->Xv();
 882             dyv =  ((AliAODMCParticle*)mcparticles0->At(dindex))->Yv();
 883           } else {
 884             dxv =  ((AliAODMCParticle*)mcparticles1->At(dindex))->Xv();
 885             dyv =  ((AliAODMCParticle*)mcparticles1->At(dindex))->Yv();
 886           }
 887           Double_t dradius = TMath::Sqrt(dxv*dxv+dyv*dyv);
 888           //printf("\tDaughter production radius = %2.2f\n",dradius);
 889           //if you convert/decay within the trackable zone, discard
 890           //from ntuple ?
 891           //      if(dradius < 200.) continue;
 892         }
 893         */
 894
 895         //find out what the ancestry of this electron is
 896         Int_t mctag = -1;
 897         Int_t input = 0;
 898         Int_t ival = ipart;
 899         if(ipart > npart0) { ival -= npart0; input = 1;}
 900         mctag = GetMCAnalysisUtils()->CheckOrigin(ival,GetReader(),input);
 901
 902         //fill ntuple
 903         if(fWriteNtuple) {
 904           fMCEleNtuple->Fill(mctag,aodprimary->Pt(),aodprimary->Phi(),aodprimary->Eta(),aodprimary->Xv(),aodprimary->Yv(),aodprimary->Zv());
 905         }
 906
 907       }
 908     }
 909   } //pure MC kine histos
 910
 911 }
 912
 913 //__________________________________________________________________
 914 Int_t AliAnaElectron::GetBtag(AliAODTrack * tr )
 915 {
 916   //This method uses the Displaced Vertex between electron-hadron
 917   //pairs and the primary vertex to determine whether an electron is
 918   //likely from a B hadron.
 919
 920   Int_t ncls1 = 0;
 921   for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;
 922   if (ncls1 < fITSCut) return 0;
 923
 924   Double_t x[3];
 925   //Note: 02-Sep-2009, Must use GetPosition, not XYZAtDCA
 926   //Bool_t gotit = tr->XYZAtDCA(x);
 927   Bool_t isNotDCA = tr->GetPosition(x);
 928   if(isNotDCA) { printf("##Problem getting impact parameter!\n"); return 0; }
 929
 930   Double_t d1 = TMath::Sqrt(x[0]*x[0] + x[1]*x[1]);
 931   if (TMath::Abs(d1)   > fImpactCut ) return 0;
 932   if (TMath::Abs(x[2]) > fImpactCut ) return 0;
 933   //printf("----- impact parameter: x=%f, y=%f, z=%f -------\n",x[0],x[1], x[2]);
 934
 935   Int_t nvtx1 = 0;
 936   Int_t nvtx2 = 0;
 937   Int_t nvtx3 = 0;
 938
 939   for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {
 940     //loop over assoc
 941     AliAODTrack* track2 = (AliAODTrack*) (GetAODCTS()->At(k2));
 942     Int_t id1 = tr->GetID();
 943     Int_t id2 = track2->GetID();
 944     if(id1 == id2) continue;
 945
 946     Int_t ncls2 = 0;
 947     for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
 948     if (ncls2 < fITSCut) continue;
 949
 950     if(track2->Pt() < fAssocPtCut) continue;
 951
 952     Double_t dphi = tr->Phi() - track2->Phi();
 953     if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
 954     if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
 955     Double_t deta = tr->Eta() - track2->Eta();
 956     Double_t dr = sqrt(deta*deta + dphi*dphi);
 957
 958     if(dr > fDrCut) continue;
 959
 960     Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);
 961     if (sDca1 > fSdcaCut) nvtx1++;
 962     Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);
 963     if (sDca2 > fSdcaCut) nvtx2++;
 964     Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);
 965     if (sDca3 > fSdcaCut) nvtx3++;
 966
 967   } //loop over hadrons
 968
 969   if(GetDebug() > 0) {
 970     if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);
 971     if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);
 972     if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);
 973   }
 974
 975   //fill QA histograms
 976   fhBtagCut1->Fill(nvtx1,tr->Pt());
 977   fhBtagCut2->Fill(nvtx2,tr->Pt());
 978   fhBtagCut3->Fill(nvtx3,tr->Pt());
 979
 980   return nvtx2;
 981 }
 982
 983 //__________________________________________________________________
 984 Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)
 985 {
 986   //Compute the signed dca between two tracks
 987   //and return the result
 988
 989   Double_t signDca=-999.;
 990   if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);
 991
 992   //=====Now calculate DCA between both tracks=======
 993   Double_t massE = 0.000511;
 994   Double_t massK = 0.493677;
 995
 996   Double_t bfield = 5.; //kG
 997   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
 998
 999   Double_t vertex[3] = {-999.,-999.,-999}; //vertex
1000   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
1001     GetReader()->GetVertex(vertex); //If only one file, get the vertex from there
1002     //FIXME:  Add a check for whether file 2 is PYTHIA or HIJING
1003     //If PYTHIA, then set the vertex from file 2, if not, use the
1004     //vertex from file 1
1005     if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);
1006   }
1007
1008   TVector3 primV(vertex[0],vertex[1],vertex[2]) ;
1009
1010   if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;
1011
1012   AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);
1013   AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);
1014
1015   Double_t xplane1 = 0.; Double_t xplane2 = 0.;
1016   Double_t pairdca = param1->GetDCA(param2,bfield,xplane1,xplane2);
1017
1018   Int_t id1 = 0, id2 = 0;
1019   AliESDv0 bvertex(*param1,id1,*param2,id2);
1020   Double_t vx,vy,vz;
1021   bvertex.GetXYZ(vx,vy,vz);
1022
1023   Double_t emom[3];
1024   Double_t hmom[3];
1025   param1->PxPyPz(emom);
1026   param2->PxPyPz(hmom);
1027   TVector3 emomAtB(emom[0],emom[1],emom[2]);
1028   TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);
1029   TVector3 secvtxpt(vx,vy,vz);
1030   TVector3 decayvector(0,0,0);
1031   decayvector = secvtxpt - primV; //decay vector from PrimVtx
1032   Double_t decaylength = decayvector.Mag();
1033
1034   if(GetDebug() > 0) {
1035     printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );
1036     printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );
1037   }
1038
1039   if (masscut<1.1) fhBtagQA1->Fill(pairdca,decaylength);
1040
1041   if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {
1042     TVector3 sumMom = emomAtB+hmomAtB;
1043     Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);
1044     Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);
1045     Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);
1046     Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
1047     Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));
1048     Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();
1049
1050     if (masscut<1.1) fhBtagQA2->Fill(sDca, mass);
1051
1052     if (mass > masscut && massPhot > 0.1) signDca = sDca;
1053
1054     if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);
1055     if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);
1056   }
1057
1058   //clean up
1059   delete param1;
1060   delete param2;
1061
1062   return signDca;
1063 }
1064
1065 //__________________________________________________________________
1066 Bool_t AliAnaElectron::IsItPhotonic(const AliAODPWG4Particle* part)
1067 {
1068   //This method checks the opening angle and invariant mass of
1069   //electron pairs to see if they are likely to be photonic electrons
1070
1071   Bool_t itIS = kFALSE;
1072
1073   Double_t massE = 0.000511;
1074   Double_t bfield = 5.; //kG
1075   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
1076
1077   Int_t pdg1 = part->GetPdg();
1078   Int_t trackId = part->GetTrackLabel(0);
1079   AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);
1080   if(!track) {
1081     if(GetDebug() > 0) printf("AliAnaElectron::IsItPhotonic - can't get the AOD Track from the particle!  Skipping the photonic check");
1082     return kFALSE; //Don't proceed because we can't get the track
1083   }
1084
1085   AliExternalTrackParam *param1 = new AliExternalTrackParam(track);
1086
1087   //Loop on stored AOD electrons and compute the angle differences and Minv
1088   for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {
1089     AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);
1090     Int_t track2Id = part2->GetTrackLabel(0);
1091     if(trackId == track2Id) continue;
1092     Int_t pdg2 = part2->GetPdg();
1093     if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;
1094     if(part2->GetDetector() != fCalorimeter) continue;
1095
1096     //JLK: Check opp. sign pairs only
1097     if(pdg1*pdg2 > 0) continue; //skip same-sign pairs
1098
1099     //propagate to common vertex and check opening angle
1100     AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(track2Id);
1101     if(!track2) {
1102       if(GetDebug() >0) printf("AliAnaElectron::IsItPhotonic - problem getting the partner track.  Continuing on to the next one");
1103       continue;
1104     }
1105     AliExternalTrackParam *param2 = new AliExternalTrackParam(track2);
1106     Int_t id1 = 0, id2 = 0;
1107     AliESDv0 photonVtx(*param1,id1,*param2,id2);
1108     Double_t vx,vy,vz;
1109     photonVtx.GetXYZ(vx,vy,vz);
1110
1111     Double_t p1mom[3];
1112     Double_t p2mom[3];
1113     param1->PxPyPz(p1mom);
1114     param2->PxPyPz(p2mom);
1115
1116     TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);
1117     TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);
1118     TVector3 sumMom = p1momAtB+p2momAtB;
1119
1120     Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);
1121     Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);
1122     Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
1123
1124     Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);
1125     fh1OpeningAngle->Fill(dphi);
1126     fh1MinvPhoton->Fill(mass);
1127
1128     if(mass < 0.1) {
1129       if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");
1130       itIS = kTRUE;
1131     }
1132
1133     //clean up
1134     delete param2;
1135
1136   }
1137
1138   delete param1;
1139   return itIS;
1140
1141 }
1142
1143 //__________________________________________________________________
1144 void AliAnaElectron::Print(const Option_t * opt) const
1145 {
1146   //Print some relevant parameters set for the analysis
1147
1148   if(! opt)
1149     return;
1150
1151   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1152   AliAnaPartCorrBaseClass::Print(" ");
1153
1154   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1155   printf("pOverE range           =     %f - %f\n",fpOverEmin,fpOverEmax);
1156   printf("residual cut           =     %f\n",fResidualCut);
1157   printf("---Btagging\n");
1158   printf("max IP-cut (e,h)       =     %f\n",fImpactCut);
1159   printf("min ITS-hits           =     %d\n",fITSCut);
1160   printf("max dR (e,h)           =     %f\n",fDrCut);
1161   printf("max pairDCA            =     %f\n",fPairDcaCut);
1162   printf("max decaylength        =     %f\n",fDecayLenCut);
1163   printf("min Associated Pt      =     %f\n",fAssocPtCut);
1164   printf("    \n") ;
1165
1166 }
1167
1168 //________________________________________________________________________
1169 void AliAnaElectron::ReadHistograms(TList* outputList)
1170 {
1171   // Needed when Terminate is executed in distributed environment
1172   // Refill analysis histograms of this class with corresponding
1173   // histograms in output list.
1174
1175   // Histograms of this analsys are kept in the same list as other
1176   // analysis, recover the position of
1177   // the first one and then add the next
1178   Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));
1179
1180   //Read histograms, must be in the same order as in
1181   //GetCreateOutputObject.
1182   fh1pOverE     = (TH1F *) outputList->At(index);
1183   fh1dR         = (TH1F *) outputList->At(index++);
1184   fh2EledEdx    = (TH2F *) outputList->At(index++);
1185   fh2MatchdEdx  = (TH2F *) outputList->At(index++);
1186
1187 }
1188
1189 //__________________________________________________________________
1190 void  AliAnaElectron::Terminate(TList* outputList)
1191 {
1192
1193   //Do some plots to end
1194   //Recover histograms from output histograms list, needed for
1195   //distributed analysis.
1196   //ReadHistograms(outputList);
1197
1198   printf(" AliAnaElectron::Terminate()  *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;
1199
1200 }
1201