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