PWG4/PartCorrDep/AliAnaElectron.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes hereby granted      *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 /* $Id: $ */\r
  16 \r
  17 //_________________________________________________________________________\r
  18 //\r
  19 // Class for the electron identification.\r
  20 // Clusters from EMCAL matched to tracks\r
  21 // and kept in the AOD. Few histograms produced.\r
  22 //\r
  23 // -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)\r
  24 //////////////////////////////////////////////////////////////////////////////\r
  25   \r
  26 // --- ROOT system --- \r
  27 #include <TH2F.h>\r
  28 #include <TH3F.h>\r
  29 #include <TParticle.h>\r
  30 #include <TNtuple.h>\r
  31 #include <TClonesArray.h>\r
  32 #include <TObjString.h>\r
  33 //#include <Riostream.h>\r
  34 \r
  35 // --- Analysis system --- \r
  36 #include "AliAnaElectron.h" \r
  37 #include "AliCaloTrackReader.h"\r
  38 #include "AliMCAnalysisUtils.h"\r
  39 #include "AliAODCaloCluster.h"\r
  40 #include "AliFidutialCut.h"\r
  41 #include "AliAODTrack.h"\r
  42 #include "AliAODPid.h"\r
  43 #include "AliCaloPID.h"\r
  44 #include "AliAODMCParticle.h"\r
  45 #include "AliStack.h"\r
  46 #include "AliExternalTrackParam.h"\r
  47 #include "AliESDv0.h"\r
  48 #include "AliESDtrack.h"\r
  49 #include "AliAODJet.h"\r
  50 #include "AliAODEvent.h"\r
  51 #include "AliGenPythiaEventHeader.h"\r
  52 \r
  53 ClassImp(AliAnaElectron)\r
  54   \r
  55 //____________________________________________________________________________\r
  56 AliAnaElectron::AliAnaElectron() \r
  57 : AliAnaPartCorrBaseClass(),fCalorimeter(""),\r
  58   fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.),\r
  59   fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),\r
  60   fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),\r
  61   fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),\r
  62   fWriteNtuple(kFALSE),\r
  63   //event QA histos\r
  64   fhImpactXY(0),fhRefMult(0),fhRefMult2(0),\r
  65   //matching checks\r
  66   fh3pOverE(0),fh3EOverp(0),fh3pOverE2(0),fh3EOverp2(0),fh3pOverE3(0),fh3EOverp3(0),\r
  67   fh2pOverE(0),fh2EOverp(0),fh2pOverE2(0),fh2EOverp2(0),\r
  68   fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),\r
  69   fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),fh2TrackPVsClusterE(0),\r
  70   fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),\r
  71   //Photonic electron checks\r
  72   fh1OpeningAngle(0),fh1MinvPhoton(0),\r
  73   //Reconstructed electrons\r
  74   fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),\r
  75   fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),\r
  76   fhPtPE(0),fhPhiPE(0),fhEtaPE(0),\r
  77   //DVM B-tagging\r
  78   fhDVMBtagCut1(0),fhDVMBtagCut2(0),fhDVMBtagCut3(0),fhDVMBtagQA1(0),fhDVMBtagQA2(0),\r
  79   fhDVMBtagQA3(0),fhDVMBtagQA4(0),fhDVMBtagQA5(0),\r
  80   //IPSig B-tagging\r
  81   fhIPSigBtagQA1(0),fhIPSigBtagQA2(0),fhTagJetPt1x4(0),fhTagJetPt2x3(0),fhTagJetPt3x2(0),\r
  82   fhePlusTagJetPt1x4(0),fhePlusTagJetPt2x3(0),fhePlusTagJetPt3x2(0),\r
  83   //B-Jet histograms\r
  84   fhJetType(0),fhLeadJetType(0),fhBJetXsiFF(0),fhBJetPtFF(0),fhBJetEtaPhi(0),\r
  85   fhNonBJetXsiFF(0),fhNonBJetPtFF(0),fhNonBJetEtaPhi(0),\r
  86   /////////////////////////////////////////////////////////////\r
  87   //Histograms that rely on MC info (not filled for real data)\r
  88   fEleNtuple(0),\r
  89   //reco electrons from various sources\r
  90   fhPhiConversion(0),fhEtaConversion(0),\r
  91   //for comparisons with tracking detectors\r
  92   fhPtTrack(0),fhPtHadron(0),fhPtNPEleTPC(0),fhPtNPEleTPCTRD(0),fhPtNPEleTTE(0),\r
  93   fhPtNPEleEMCAL(0),fhPtNPEBHadron(0),\r
  94   //for computing efficiency of B-jet tags\r
  95   fhBJetPt1x4(0),fhBJetPt2x3(0),fhBJetPt3x2(0),\r
  96   fhFakeJetPt1x4(0),fhFakeJetPt2x3(0),fhFakeJetPt3x2(0),fhDVMJet(0),\r
  97   //MC rate histograms/ntuple\r
  98   fMCEleNtuple(0),fhMCBJetElePt(0),fhMCBHadronElePt(0),fhPtMCHadron(0),fhPtMCElectron(0),\r
  99   fhMCXYConversion(0),fhMCRadPtConversion(0)\r
 100 {\r
 101   //default ctor\r
 102   \r
 103   //Initialize parameters\r
 104   InitParameters();\r
 105 \r
 106 }\r
 107 \r
 108 //____________________________________________________________________________\r
 109 AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) \r
 110   : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),\r
 111     fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),\r
 112     fResidualCut(g.fResidualCut),fMinClusEne(g.fMinClusEne),\r
 113     fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),\r
 114     fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),\r
 115     fNTagTrkCut(g.fNTagTrkCut),fIPSigCut(g.fIPSigCut),\r
 116     fJetEtaCut(g.fJetEtaCut),fJetPhiMin(g.fJetPhiMin),fJetPhiMax(g.fJetPhiMax),\r
 117     fWriteNtuple(g.fWriteNtuple),\r
 118     //event QA histos\r
 119     fhImpactXY(g.fhImpactXY),fhRefMult(g.fhRefMult),fhRefMult2(g.fhRefMult2),\r
 120     //matching checks\r
 121     fh3pOverE(g.fh3pOverE),fh3EOverp(g.fh3EOverp),\r
 122     fh3pOverE2(g.fh3pOverE2),fh3EOverp2(g.fh3EOverp2),\r
 123     fh3pOverE3(g.fh3pOverE3),fh3EOverp3(g.fh3EOverp3),\r
 124     fh2pOverE(g.fh2pOverE),fh2EOverp(g.fh2EOverp),\r
 125     fh2pOverE2(g.fh2pOverE2),fh2EOverp2(g.fh2EOverp2),\r
 126     fh1dR(g.fh1dR),fh2EledEdx(g.fh2EledEdx),\r
 127     fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),\r
 128     fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),\r
 129     fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),\r
 130     fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),\r
 131     //Photonic electron checks\r
 132     fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),\r
 133     //Reconstructed electrons\r
 134     fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron),\r
 135     fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE),\r
 136     fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE),\r
 137     //DVM B-tagging\r
 138     fhDVMBtagCut1(g.fhDVMBtagCut1),fhDVMBtagCut2(g.fhDVMBtagCut2),fhDVMBtagCut3(g.fhDVMBtagCut3),\r
 139     fhDVMBtagQA1(g.fhDVMBtagQA1),fhDVMBtagQA2(g.fhDVMBtagQA2),\r
 140     fhDVMBtagQA3(g.fhDVMBtagQA3),fhDVMBtagQA4(g.fhDVMBtagQA4),fhDVMBtagQA5(g.fhDVMBtagQA5),\r
 141     //IPSig B-tagging\r
 142     fhIPSigBtagQA1(g.fhIPSigBtagQA1),fhIPSigBtagQA2(g.fhIPSigBtagQA2),\r
 143     fhTagJetPt1x4(g.fhTagJetPt1x4),fhTagJetPt2x3(g.fhTagJetPt2x3),fhTagJetPt3x2(g.fhTagJetPt3x2),\r
 144     fhePlusTagJetPt1x4(g.fhePlusTagJetPt1x4),fhePlusTagJetPt2x3(g.fhePlusTagJetPt2x3),\r
 145     fhePlusTagJetPt3x2(g.fhePlusTagJetPt3x2),\r
 146     //B-Jet histograms\r
 147     fhJetType(g.fhJetType),fhLeadJetType(g.fhLeadJetType),fhBJetXsiFF(g.fhBJetXsiFF),\r
 148     fhBJetPtFF(g.fhBJetPtFF),fhBJetEtaPhi(g.fhBJetEtaPhi),fhNonBJetXsiFF(g.fhNonBJetXsiFF),\r
 149     fhNonBJetPtFF(g.fhNonBJetPtFF),fhNonBJetEtaPhi(g.fhNonBJetEtaPhi),\r
 150     /////////////////////////////////////////////////////////////\r
 151     //Histograms that rely on MC info (not filled for real data)\r
 152     fEleNtuple(g.fEleNtuple),\r
 153     //reco electrons from various sources\r
 154     fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),\r
 155     //for comparisons with tracking detectors\r
 156     fhPtTrack(g.fhPtTrack),fhPtHadron(g.fhPtHadron),fhPtNPEleTPC(g.fhPtNPEleTPC),\r
 157     fhPtNPEleTPCTRD(g.fhPtNPEleTPCTRD),fhPtNPEleTTE(g.fhPtNPEleTTE),\r
 158     fhPtNPEleEMCAL(g.fhPtNPEleEMCAL),fhPtNPEBHadron(g.fhPtNPEBHadron),\r
 159     //for computing efficiency of B-jet tags\r
 160     fhBJetPt1x4(g.fhBJetPt1x4),fhBJetPt2x3(g.fhBJetPt2x3),\r
 161     fhBJetPt3x2(g.fhBJetPt3x2),\r
 162     fhFakeJetPt1x4(g.fhFakeJetPt1x4),fhFakeJetPt2x3(g.fhBJetPt2x3),\r
 163     fhFakeJetPt3x2(g.fhFakeJetPt3x2),fhDVMJet(g.fhDVMJet),\r
 164     //MC rate histograms/ntuple\r
 165     fMCEleNtuple(g.fMCEleNtuple),fhMCBJetElePt(g.fhMCBJetElePt),\r
 166     fhMCBHadronElePt(g.fhMCBHadronElePt),\r
 167     fhPtMCHadron(g.fhPtMCHadron),fhPtMCElectron(g.fhPtMCElectron),\r
 168     fhMCXYConversion(g.fhMCXYConversion),fhMCRadPtConversion(g.fhMCRadPtConversion)\r
 169 {\r
 170   // cpy ctor\r
 171   \r
 172 }\r
 173 \r
 174 //_________________________________________________________________________\r
 175 AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)\r
 176 {\r
 177   // assignment operator\r
 178   \r
 179   if(&g == this) return *this;\r
 180   fCalorimeter = g.fCalorimeter;\r
 181   fpOverEmin = g.fpOverEmin;\r
 182   fpOverEmax = g.fpOverEmax;\r
 183   fResidualCut = g.fResidualCut;\r
 184   fMinClusEne = g.fMinClusEne;\r
 185   fDrCut = g.fDrCut;\r
 186   fPairDcaCut = g.fPairDcaCut;\r
 187   fDecayLenCut = g.fDecayLenCut;\r
 188   fImpactCut = g.fImpactCut;\r
 189   fAssocPtCut = g.fAssocPtCut;\r
 190   fMassCut = g.fMassCut;\r
 191   fSdcaCut = g.fSdcaCut;\r
 192   fITSCut = g.fITSCut;\r
 193   fNTagTrkCut = g.fNTagTrkCut;\r
 194   fIPSigCut = g.fIPSigCut;\r
 195   fJetEtaCut = g.fJetEtaCut;\r
 196   fJetPhiMin = g.fJetPhiMin;\r
 197   fJetPhiMax = g.fJetPhiMax;\r
 198   fWriteNtuple = g.fWriteNtuple;\r
 199   //event QA histos\r
 200   fhImpactXY = g.fhImpactXY;\r
 201   fhRefMult  = g.fhRefMult;\r
 202   fhRefMult2 = g.fhRefMult2;\r
 203   //matching checks\r
 204   fh3pOverE  = g.fh3pOverE;\r
 205   fh3EOverp  = g.fh3EOverp;\r
 206   fh3pOverE2 = g.fh3pOverE2;\r
 207   fh3EOverp2 = g.fh3EOverp2;\r
 208   fh3pOverE3 = g.fh3pOverE3;\r
 209   fh3EOverp3 = g.fh3EOverp3;\r
 210   fh2pOverE  = g.fh2pOverE;\r
 211   fh2EOverp  = g.fh2EOverp;\r
 212   fh2pOverE2 = g.fh2pOverE2;\r
 213   fh2EOverp2 = g.fh2EOverp2;\r
 214   fh1dR     = g.fh1dR;\r
 215   fh2EledEdx = g.fh2EledEdx;\r
 216   fh2MatchdEdx = g.fh2MatchdEdx;\r
 217   fh2dEtadPhi = g.fh2dEtadPhi;\r
 218   fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;\r
 219   fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;\r
 220   fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;\r
 221   fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;\r
 222   fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;\r
 223   fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;\r
 224   //Photonic electron checks\r
 225   fh1OpeningAngle = g.fh1OpeningAngle;\r
 226   fh1MinvPhoton = g.fh1MinvPhoton;\r
 227   //Reconstructed electrons\r
 228   fhPtElectron = g.fhPtElectron; \r
 229   fhPhiElectron = g.fhPhiElectron; \r
 230   fhEtaElectron = g.fhEtaElectron; \r
 231   fhPtNPE = g.fhPtNPE;\r
 232   fhPhiNPE = g.fhPhiNPE;\r
 233   fhEtaNPE = g.fhEtaNPE; \r
 234   fhPtPE = g.fhPtPE;\r
 235   fhPhiPE = g.fhPhiPE;\r
 236   fhEtaPE = g.fhEtaPE; \r
 237   //DVM B-tagging\r
 238   fhDVMBtagCut1 = g.fhDVMBtagCut1;\r
 239   fhDVMBtagCut2 = g.fhDVMBtagCut2; \r
 240   fhDVMBtagCut3 = g.fhDVMBtagCut3; \r
 241   fhDVMBtagQA1 = g.fhDVMBtagQA1; \r
 242   fhDVMBtagQA2 = g.fhDVMBtagQA2; \r
 243   fhDVMBtagQA3 = g.fhDVMBtagQA3; \r
 244   fhDVMBtagQA4 = g.fhDVMBtagQA4; \r
 245   fhDVMBtagQA5 = g.fhDVMBtagQA5; \r
 246   //IPSig B-tagging\r
 247   fhIPSigBtagQA1 = g.fhIPSigBtagQA1; \r
 248   fhIPSigBtagQA2 = g.fhIPSigBtagQA2; \r
 249   fhTagJetPt1x4 = g.fhTagJetPt1x4; \r
 250   fhTagJetPt2x3 = g.fhTagJetPt2x3; \r
 251   fhTagJetPt3x2 = g.fhTagJetPt3x2; \r
 252   fhePlusTagJetPt1x4 = g.fhePlusTagJetPt1x4; \r
 253   fhePlusTagJetPt2x3 = g.fhePlusTagJetPt2x3; \r
 254   fhePlusTagJetPt3x2 = g.fhePlusTagJetPt3x2; \r
 255   //B-Jet histograms\r
 256   fhJetType = g.fhJetType; \r
 257   fhLeadJetType = g.fhLeadJetType; \r
 258   fhBJetXsiFF = g.fhBJetXsiFF; \r
 259   fhBJetPtFF = g.fhBJetPtFF; \r
 260   fhBJetEtaPhi = g.fhBJetEtaPhi; \r
 261   fhNonBJetXsiFF = g.fhNonBJetXsiFF; \r
 262   fhNonBJetPtFF = g.fhNonBJetPtFF; \r
 263   fhNonBJetEtaPhi = g.fhNonBJetEtaPhi; \r
 264   /////////////////////////////////////////////////////////////\r
 265   //Histograms that rely on MC info (not filled for real data)\r
 266   fEleNtuple = g.fEleNtuple; \r
 267   //reco electrons from various sources\r
 268   fhPhiConversion = g.fhPhiConversion; \r
 269   fhEtaConversion = g.fhEtaConversion;\r
 270   //for comparisons with tracking detectors\r
 271   fhPtTrack = g.fhPtTrack;\r
 272   fhPtHadron = g.fhPtHadron; fhPtNPEleTPC = g.fhPtNPEleTPC; \r
 273   fhPtNPEleTPCTRD = g.fhPtNPEleTPCTRD; fhPtNPEleTTE = g.fhPtNPEleTTE; \r
 274   fhPtNPEleEMCAL = g.fhPtNPEleEMCAL; fhPtNPEBHadron = g.fhPtNPEBHadron;\r
 275   //for computing efficiency of B-jet tags\r
 276   fhBJetPt1x4 = g.fhBJetPt1x4; fhBJetPt2x3 = g.fhBJetPt2x3; \r
 277   fhBJetPt3x2 = g.fhBJetPt3x2;\r
 278   fhFakeJetPt1x4 = g.fhFakeJetPt1x4; fhFakeJetPt2x3 = g.fhFakeJetPt2x3; \r
 279   fhFakeJetPt3x2 = g.fhFakeJetPt3x2; fhDVMJet = g.fhDVMJet;\r
 280   //MC rate histograms/ntuple\r
 281   fMCEleNtuple = g.fMCEleNtuple; fhMCBJetElePt = g.fhMCBJetElePt; \r
 282   fhMCBHadronElePt = g.fhMCBHadronElePt;\r
 283   fhPtMCHadron = g.fhPtMCHadron; fhPtMCElectron = g.fhPtMCElectron; \r
 284   fhMCXYConversion = g.fhMCXYConversion;\r
 285   fhMCRadPtConversion = g.fhMCRadPtConversion;\r
 286 \r
 287   return *this;\r
 288   \r
 289 }\r
 290 \r
 291 //____________________________________________________________________________\r
 292 AliAnaElectron::~AliAnaElectron() \r
 293 {\r
 294   //dtor\r
 295 \r
 296 }\r
 297 \r
 298 \r
 299 //________________________________________________________________________\r
 300 TList *  AliAnaElectron::GetCreateOutputObjects()\r
 301 {  \r
 302   // Create histograms to be saved in output file and \r
 303   // store them in outputContainer\r
 304   TList * outputContainer = new TList() ; \r
 305   outputContainer->SetName("ElectronHistos") ; \r
 306 \r
 307   Int_t nptbins  = GetHistoNPtBins();\r
 308   Int_t nphibins = GetHistoNPhiBins();\r
 309   Int_t netabins = GetHistoNEtaBins();\r
 310   Float_t ptmax  = GetHistoPtMax();\r
 311   Float_t phimax = GetHistoPhiMax();\r
 312   Float_t etamax = GetHistoEtaMax();\r
 313   Float_t ptmin  = GetHistoPtMin();\r
 314   Float_t phimin = GetHistoPhiMin();\r
 315   Float_t etamin = GetHistoEtaMin();    \r
 316 \r
 317   //event QA\r
 318   fhImpactXY = new TH1F("hImpactXY","Impact parameter for all tracks",200,-10,10.);\r
 319   fhRefMult = new TH1F("hRefMult" ,"refmult QA: " ,5000,0,5000);\r
 320   fhRefMult2  = new TH1F("hRefMult2" ,"refmult2 QA: " ,5000,0,5000);\r
 321 \r
 322   outputContainer->Add(fhImpactXY);\r
 323   outputContainer->Add(fhRefMult);\r
 324   outputContainer->Add(fhRefMult2);\r
 325   \r
 326   //matching checks\r
 327   fh3pOverE  = new TH3F("h3pOverE"    ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
 328   fh3EOverp  = new TH3F("h3EOverp"    ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
 329   fh3pOverE2 = new TH3F("h3pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
 330   fh3EOverp2 = new TH3F("h3EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
 331   fh3pOverE3 = new TH3F("h3pOverE_Tpc","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
 332   fh3EOverp3 = new TH3F("h3EOverp_Tpc","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
 333   fh2pOverE  = new TH2F("h2pOverE"    ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);\r
 334   fh2EOverp  = new TH2F("h2EOverp"    ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );\r
 335   fh2pOverE2 = new TH2F("h2pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);\r
 336   fh2EOverp2 = new TH2F("h2EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );\r
 337 \r
 338   fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi());\r
 339   fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.);\r
 340   fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.);\r
 341   fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
 342   fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
 343   fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
 344 \r
 345   fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
 346   fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
 347   fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax);\r
 348   fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax);\r
 349 \r
 350   outputContainer->Add(fh3pOverE) ;\r
 351   outputContainer->Add(fh3EOverp) ;\r
 352   outputContainer->Add(fh3pOverE2) ;\r
 353   outputContainer->Add(fh3EOverp2) ;\r
 354   outputContainer->Add(fh3pOverE3) ;\r
 355   outputContainer->Add(fh3EOverp3) ;\r
 356   outputContainer->Add(fh2pOverE) ;\r
 357   outputContainer->Add(fh2EOverp) ;\r
 358   outputContainer->Add(fh2pOverE2) ;\r
 359   outputContainer->Add(fh2EOverp2) ;\r
 360   outputContainer->Add(fh1dR) ; \r
 361   outputContainer->Add(fh2EledEdx) ;\r
 362   outputContainer->Add(fh2MatchdEdx) ;\r
 363   outputContainer->Add(fh2dEtadPhi) ;\r
 364   outputContainer->Add(fh2dEtadPhiMatched) ;\r
 365   outputContainer->Add(fh2dEtadPhiUnmatched) ;\r
 366   outputContainer->Add(fh2TrackPVsClusterE) ;\r
 367   outputContainer->Add(fh2TrackPtVsClusterE) ;\r
 368   outputContainer->Add(fh2TrackPhiVsClusterPhi) ;\r
 369   outputContainer->Add(fh2TrackEtaVsClusterEta) ;\r
 370   \r
 371   //photonic electron checks\r
 372   fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between e+e- pairs",100,0.,TMath::Pi());\r
 373   fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of e+e- pairs",200,0.,2.);\r
 374 \r
 375   outputContainer->Add(fh1OpeningAngle);\r
 376   outputContainer->Add(fh1MinvPhoton);\r
 377 \r
 378   //Reconstructed electrons\r
 379   fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);\r
 380   fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
 381   fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
 382   fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);\r
 383   fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
 384   fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
 385   fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);\r
 386   fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
 387   fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
 388 \r
 389   outputContainer->Add(fhPtElectron) ; \r
 390   outputContainer->Add(fhPhiElectron) ; \r
 391   outputContainer->Add(fhEtaElectron) ;\r
 392   outputContainer->Add(fhPtNPE) ; \r
 393   outputContainer->Add(fhPhiNPE) ; \r
 394   outputContainer->Add(fhEtaNPE) ;\r
 395   outputContainer->Add(fhPtPE) ; \r
 396   outputContainer->Add(fhPhiPE) ; \r
 397   outputContainer->Add(fhEtaPE) ;\r
 398 \r
 399   //B-tagging\r
 400   fhDVMBtagCut1 = new TH2F("hdvmbtag_cut1","DVM B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);\r
 401   fhDVMBtagCut2 = new TH2F("hdvmbtag_cut2","DVM B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);\r
 402   fhDVMBtagCut3 = new TH2F("hdvmbtag_cut3","DVM B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);\r
 403   fhDVMBtagQA1  = new TH2F("hdvmbtag_qa1" ,"DVM B-tag QA: pairDCA vs length", 100,0,0.2 ,100,0,1.0);\r
 404   fhDVMBtagQA2  = new TH2F("hdvmbtag_qa2" ,"DVM B-tag QA: signDCA vs mass"  , 200,-0.5,0.5 ,100,0,10);\r
 405   fhDVMBtagQA3  = new TH1F("hdvmbtag_qa3" ,"DVM B-tag QA: ITS-Hits electron" ,7,0,7);\r
 406   fhDVMBtagQA4  = new TH1F("hdvmbtag_qa4" ,"DVM B-tag QA: IP d electron" ,200,-3,3);\r
 407   fhDVMBtagQA5  = new TH1F("hdvmbtag_qa5" ,"DVM B-tag QA: IP z electron" ,200,-3,3);\r
 408 \r
 409   outputContainer->Add(fhDVMBtagCut1) ;\r
 410   outputContainer->Add(fhDVMBtagCut2) ;\r
 411   outputContainer->Add(fhDVMBtagCut3) ;\r
 412   outputContainer->Add(fhDVMBtagQA1) ;\r
 413   outputContainer->Add(fhDVMBtagQA2) ;\r
 414   outputContainer->Add(fhDVMBtagQA3) ;\r
 415   outputContainer->Add(fhDVMBtagQA4) ;\r
 416   outputContainer->Add(fhDVMBtagQA5) ;\r
 417 \r
 418   //IPSig B-tagging\r
 419   fhIPSigBtagQA1  = new TH1F("hipsigbtag_qa1" ,"IPSig B-tag QA: # tag tracks", 20,0,20);\r
 420   fhIPSigBtagQA2  = new TH1F("hipsigbtag_qa2" ,"IPSig B-tag QA: IP significance", 200,-10.,10.);\r
 421   fhTagJetPt1x4 = new TH1F("hTagJetPt1x4","tagged jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
 422   fhTagJetPt2x3 = new TH1F("hTagJetPt2x3","tagged jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
 423   fhTagJetPt3x2 = new TH1F("hTagJetPt3x2","tagged jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
 424   fhePlusTagJetPt1x4 = new TH1F("hePlusTagJetPt1x4","tagged eJet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
 425   fhePlusTagJetPt2x3 = new TH1F("hePlusTagJetPt2x3","tagged eJet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
 426   fhePlusTagJetPt3x2 = new TH1F("hePlusTagJetPt3x2","tagged eJet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
 427 \r
 428   outputContainer->Add(fhIPSigBtagQA1) ;\r
 429   outputContainer->Add(fhIPSigBtagQA2) ;\r
 430   outputContainer->Add(fhTagJetPt1x4);\r
 431   outputContainer->Add(fhTagJetPt2x3);\r
 432   outputContainer->Add(fhTagJetPt3x2);\r
 433   outputContainer->Add(fhePlusTagJetPt1x4);\r
 434   outputContainer->Add(fhePlusTagJetPt2x3);\r
 435   outputContainer->Add(fhePlusTagJetPt3x2);\r
 436 \r
 437   //B-Jet histograms\r
 438   fhJetType = new TH2F("hJetType","# jets passing each tag method vs jet pt",15,0,15,300,0.,300.);\r
 439   fhLeadJetType = new TH2F("hLeadJetType","# leading jets passing each tag method vs jet pt",15,0,15,300,0.,300.);\r
 440   fhBJetXsiFF = new TH2F("hBJetXsiFF","B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);\r
 441   fhBJetPtFF = new TH2F("hBJetPtFF","B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);\r
 442   fhBJetEtaPhi = new TH2F("hBJetEtaPhi","B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);\r
 443   fhNonBJetXsiFF = new TH2F("hNonBJetXsiFF","Non B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);\r
 444   fhNonBJetPtFF = new TH2F("hNonBJetPtFF","Non B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);\r
 445   fhNonBJetEtaPhi = new TH2F("hNonBJetEtaPhi","Non B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);\r
 446 \r
 447   outputContainer->Add(fhJetType);\r
 448   outputContainer->Add(fhLeadJetType);\r
 449   outputContainer->Add(fhBJetXsiFF);\r
 450   outputContainer->Add(fhBJetPtFF);\r
 451   outputContainer->Add(fhBJetEtaPhi);\r
 452   outputContainer->Add(fhNonBJetXsiFF);\r
 453   outputContainer->Add(fhNonBJetPtFF);\r
 454   outputContainer->Add(fhNonBJetEtaPhi);\r
 455 \r
 456   //Histograms that use MC information\r
 457   if(IsDataMC()){\r
 458 \r
 459     //electron ntuple for further analysis\r
 460     if(fWriteNtuple) {\r
 461       fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");\r
 462       outputContainer->Add(fEleNtuple) ;\r
 463     }\r
 464 \r
 465     //electrons from various MC sources\r
 466     fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
 467     fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. pT",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
 468 \r
 469     outputContainer->Add(fhPhiConversion);\r
 470     outputContainer->Add(fhEtaConversion);\r
 471 \r
 472     //Bins along y-axis are:  0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,\r
 473     //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown, 8 - misidentified\r
 474 \r
 475     //histograms for comparison to tracking detectors\r
 476     fhPtTrack  = new TH2F("hPtTrack","Track w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 477     fhPtHadron = new TH2F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 478     fhPtNPEleTPC = new TH2F("hPtNPEleTPC","Non-phot. Electrons identified by TPC w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 479     fhPtNPEleTPCTRD = new TH2F("hPtNPEleTPCTRD","Non-phot. Electrons identified by TPC+TRD w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 480     fhPtNPEleTTE = new TH2F("hPtNPEleTTE","Non-phot. Electrons identified by TPC+TRD+EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);    \r
 481     fhPtNPEleEMCAL = new TH2F("hPtNPEleEMCAL","Non-phot. Electrons identified by EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 482     fhPtNPEBHadron = new TH2F("hPtNPEBHadron","Non-phot. b-electrons (TPC+TRD+EMCAL) vs B-hadron pt w/in EMCAL acceptance",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
 483 \r
 484     outputContainer->Add(fhPtTrack);\r
 485     outputContainer->Add(fhPtHadron);\r
 486     outputContainer->Add(fhPtNPEleTPC);\r
 487     outputContainer->Add(fhPtNPEleTPCTRD);\r
 488     outputContainer->Add(fhPtNPEleTTE);\r
 489     outputContainer->Add(fhPtNPEleEMCAL);\r
 490     outputContainer->Add(fhPtNPEBHadron);\r
 491 \r
 492     //for computing efficiency of IPSig tag\r
 493     fhBJetPt1x4 = new TH1F("hBJetPt1x4","tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
 494     fhBJetPt2x3 = new TH1F("hBJetPt2x3","tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
 495     fhBJetPt3x2 = new TH1F("hBJetPt3x2","tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
 496     fhFakeJetPt1x4 = new TH1F("hFakeJetPt1x4","fake tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
 497     fhFakeJetPt2x3 = new TH1F("hFakeJetPt2x3","fake tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
 498     fhFakeJetPt3x2 = new TH1F("hFakeJetPt3x2","fake tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
 499     fhDVMJet = new TH2F("hDVM_algo","# DVM jets passing vs Mc-Bjet",10,0,10,300,0.,300.);\r
 500 \r
 501     outputContainer->Add(fhBJetPt1x4);\r
 502     outputContainer->Add(fhBJetPt2x3);\r
 503     outputContainer->Add(fhBJetPt3x2);\r
 504     outputContainer->Add(fhFakeJetPt1x4);\r
 505     outputContainer->Add(fhFakeJetPt2x3);\r
 506     outputContainer->Add(fhFakeJetPt3x2);\r
 507     outputContainer->Add(fhDVMJet);\r
 508 \r
 509     //MC Only histograms\r
 510     \r
 511     //MC ele ntuple for further analysis\r
 512     if(fWriteNtuple) {\r
 513       fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");\r
 514       outputContainer->Add(fMCEleNtuple) ;\r
 515     }\r
 516 \r
 517     fhMCBJetElePt = new TH2F("hMCBJetElePt","MC B-jet pT vs. electron pT",300,0.,300.,300,0.,300.);\r
 518     fhMCBHadronElePt = new TH2F("hMCBHadronElePt","MC B-hadron pT vs. electron pT",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
 519     fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);\r
 520 \r
 521     //Bins along y-axis are:  0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,\r
 522     //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown\r
 523     fhPtMCElectron = new TH2F("hPtMCElectron","MC electrons from various sources w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
 524 \r
 525     fhMCXYConversion = new TH2F("hMCXYConversion","XvsY of conversion electrons",400,-400.,400.,400,-400.,400.);\r
 526     fhMCRadPtConversion = new TH2F("hMCRadPtConversion","Radius vs pT of conversion electrons",200,0.,400.,nptbins,ptmin,ptmax);\r
 527 \r
 528     outputContainer->Add(fhMCBJetElePt);\r
 529     outputContainer->Add(fhMCBHadronElePt);\r
 530     outputContainer->Add(fhPtMCHadron);\r
 531     outputContainer->Add(fhPtMCElectron);\r
 532     outputContainer->Add(fhMCXYConversion);\r
 533     outputContainer->Add(fhMCRadPtConversion);\r
 534 \r
 535   }//Histos with MC\r
 536   \r
 537   //Save parameters used for analysis\r
 538   TString parList ; //this will be list of parameters used for this analysis.\r
 539   char onePar[500] ;\r
 540   \r
 541   sprintf(onePar,"--- AliAnaElectron ---\n") ;\r
 542   parList+=onePar ;     \r
 543   sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ;\r
 544   parList+=onePar ;  \r
 545   sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ;\r
 546   parList+=onePar ;  \r
 547   sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ;\r
 548   parList+=onePar ;  \r
 549   sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ;\r
 550   parList+=onePar ;  \r
 551   sprintf(onePar,"fMinClusEne: %f\n",fMinClusEne) ;\r
 552   parList+=onePar ;\r
 553   sprintf(onePar,"---DVM Btagging\n");\r
 554   parList+=onePar ;\r
 555   sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut);\r
 556   parList+=onePar ;\r
 557   sprintf(onePar,"min ITS-hits: %d\n",fITSCut);\r
 558   parList+=onePar ;\r
 559   sprintf(onePar,"max dR (e,h): %f\n",fDrCut);\r
 560   parList+=onePar ;\r
 561   sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut);\r
 562   parList+=onePar ;\r
 563   sprintf(onePar,"max decaylength: %f\n",fDecayLenCut);\r
 564   parList+=onePar ;\r
 565   sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut);\r
 566   parList+=onePar ;\r
 567   sprintf(onePar,"---IPSig Btagging\n");\r
 568   parList+=onePar ;\r
 569   sprintf(onePar,"min tag track: %d\n",fNTagTrkCut);\r
 570   parList+=onePar ;\r
 571   sprintf(onePar,"min IP significance: %f\n",fIPSigCut);\r
 572   parList+=onePar ;\r
 573 \r
 574   //Get parameters set in base class.\r
 575   parList += GetBaseParametersList() ;\r
 576   \r
 577   //Get parameters set in FidutialCut class (not available yet)\r
 578   //parlist += GetFidCut()->GetFidCutParametersList() \r
 579   \r
 580   TObjString *oString= new TObjString(parList) ;\r
 581   outputContainer->Add(oString);\r
 582   \r
 583   return outputContainer ;\r
 584   \r
 585 }\r
 586 \r
 587 //____________________________________________________________________________\r
 588 void AliAnaElectron::Init()\r
 589 {\r
 590 \r
 591   //do some initialization\r
 592   if(fCalorimeter == "PHOS") {\r
 593     printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");\r
 594     fCalorimeter = "EMCAL";\r
 595   }\r
 596   if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){\r
 597     printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");\r
 598     abort();\r
 599   }\r
 600 \r
 601 }\r
 602 \r
 603 \r
 604 //____________________________________________________________________________\r
 605 void AliAnaElectron::InitParameters()\r
 606 {\r
 607   \r
 608   //Initialize the parameters of the analysis.\r
 609   SetOutputAODClassName("AliAODPWG4Particle");\r
 610   SetOutputAODName("PWG4Particle");\r
 611 \r
 612   AddToHistogramsName("AnaElectron_");\r
 613 \r
 614   fCalorimeter = "EMCAL" ;\r
 615   fpOverEmin = 0.5;\r
 616   fpOverEmax = 1.2;\r
 617   fResidualCut = 0.02;\r
 618   fMinClusEne = 4.0;\r
 619   //DVM B-tagging\r
 620   fDrCut       = 1.0; \r
 621   fPairDcaCut  = 0.02;\r
 622   fDecayLenCut = 1.0;\r
 623   fImpactCut   = 0.5;\r
 624   fAssocPtCut  = 1.0;\r
 625   fMassCut     = 1.5;\r
 626   fSdcaCut     = 0.1;\r
 627   fITSCut      = 4;\r
 628   //IPSig B-tagging\r
 629   fNTagTrkCut  = 2;\r
 630   fIPSigCut    = 3.0;\r
 631   //Jet fiducial cuts\r
 632   fJetEtaCut = 0.3;\r
 633   fJetPhiMin = 1.8;\r
 634   fJetPhiMax = 2.9;\r
 635 }\r
 636 \r
 637 //__________________________________________________________________\r
 638 void  AliAnaElectron::MakeAnalysisFillAOD() \r
 639 {\r
 640   //\r
 641   // Do analysis and fill aods with electron candidates\r
 642   // These AODs will be used to do subsequent histogram filling\r
 643   //\r
 644   // Also fill some QA histograms\r
 645   //\r
 646 \r
 647   TObjArray *cl = new TObjArray();\r
 648 \r
 649   Double_t bfield = 0.;\r
 650   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
 651 \r
 652   //Select the calorimeter of the electron\r
 653   if(fCalorimeter != "EMCAL") {\r
 654     printf("This class not yet implemented for PHOS\n");\r
 655     abort();\r
 656   }\r
 657   cl = GetAODEMCAL();\r
 658   \r
 659   ////////////////////////////////////////////////\r
 660   //Start from tracks and get associated clusters \r
 661   ////////////////////////////////////////////////\r
 662   if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;\r
 663   Int_t ntracks = GetAODCTS()->GetEntriesFast();\r
 664   Int_t refmult = 0; Int_t refmult2 = 0;\r
 665   if(GetDebug() > 0)\r
 666     printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);\r
 667 \r
 668   //Unfortunately, AliAODTracks don't have associated EMCAL clusters.\r
 669   //we have to redo track-matching, I guess\r
 670   Int_t iCluster = -999;\r
 671   Int_t bt = 0; //counter for event b-tags\r
 672 \r
 673   for (Int_t itrk =  0; itrk <  ntracks; itrk++) {////////////// track loop\r
 674     iCluster = -999; //start with no match\r
 675     AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;\r
 676     if (TMath::Abs(track->Eta())< 0.5) refmult++;\r
 677     Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
 678     Bool_t dcaOkay = GetDCA(track,imp,cov);  //homegrown dca calculation until AOD is fixed\r
 679     if(!dcaOkay) printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d.  Skipping it...\n",itrk);\r
 680     if(TMath::Abs(track->Eta())< 0.5 && TMath::Abs(imp[0])<1.0 && TMath::Abs(imp[1])<1.0) refmult2++;\r
 681     fhImpactXY->Fill(imp[0]);\r
 682 \r
 683     //JLK CHECK\r
 684     //AliESDtrack esdTrack(track);\r
 685     //Double_t tpcpid[AliPID::kSPECIES];\r
 686     //esdTrack.GetTPCpid(tpcpid);\r
 687     //Double_t eProb = tpcpid[AliPID::kElectron];\r
 688     //if(eProb > 0) printf("<%d> ESD eProb = %2.2f\n",itrk,eProb);\r
 689 \r
 690     AliAODPid* pid = (AliAODPid*) track->GetDetPid();\r
 691     if(pid == 0) {\r
 692       if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);\r
 693       continue;\r
 694     } else {\r
 695       Double_t emcpos[3];\r
 696       pid->GetEMCALPosition(emcpos);\r
 697       Double_t emcmom[3];\r
 698       pid->GetEMCALMomentum(emcmom);\r
 699       \r
 700       TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);\r
 701       TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);\r
 702       Double_t tphi = pos.Phi();\r
 703       Double_t teta = pos.Eta();\r
 704       Double_t tmom = mom.Mag();\r
 705       \r
 706       TLorentzVector mom2(mom,0.);\r
 707       Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter) ;\r
 708       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);\r
 709       if(mom.Pt() > GetMinPt() && in) {\r
 710         \r
 711         Double_t dEdx = pid->GetTPCsignal();\r
 712                 \r
 713         Int_t ntot = cl->GetEntriesFast();\r
 714         Double_t res = 999.;\r
 715         Double_t pOverE = -999.;\r
 716         \r
 717         Int_t pidProb = track->GetMostProbablePID();\r
 718         Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;\r
 719         Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;\r
 720         Bool_t trkChgHad = kFALSE; if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) trkChgHad = kTRUE;\r
 721 \r
 722         Int_t tmctag = -1;\r
 723 \r
 724         //Check against V0 for conversion, only if it is flagged as electron\r
 725         Bool_t photonic = kFALSE;\r
 726         if(tpcEle || trkEle) photonic = PhotonicV0(itrk);\r
 727         if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),0); //0 = no MC info\r
 728         if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),0); //0 = no MC info\r
 729 \r
 730         if(trkChgHad) fhPtHadron->Fill(track->Pt(),0); //0 = no MC info\r
 731         if(IsDataMC()) {\r
 732           //Input from second AOD?\r
 733           Int_t input = 0;\r
 734           if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;\r
 735           tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);\r
 736 \r
 737           if(trkChgHad) fhPtHadron->Fill(track->Pt(),GetMCSource(tmctag));\r
 738           if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),GetMCSource(tmctag));\r
 739           if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),GetMCSource(tmctag));\r
 740           fhPtTrack->Fill(track->Pt(),GetMCSource(tmctag));\r
 741         }\r
 742 \r
 743         Bool_t emcEle = kFALSE;      \r
 744         //For tracks in EMCAL acceptance, pair them with all clusters\r
 745         //and fill the dEta vs dPhi for these pairs:\r
 746 \r
 747         Double_t minR  = 99;\r
 748         Double_t minPe =-1;\r
 749         Double_t minEp =-1;\r
 750         Double_t minMult = -1;\r
 751         Double_t minPt   = -1;\r
 752 \r
 753         for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
 754           AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));\r
 755           if(!clus) continue;\r
 756 \r
 757           //As of 11-Oct-2009\r
 758           //only select "good" clusters   \r
 759           if (clus->GetNCells()       < 2    ) continue;\r
 760           if (clus->GetNCells()       > 30   ) continue;\r
 761           if (clus->E()               < fMinClusEne ) continue;\r
 762           if (clus->GetDispersion()   > 1    ) continue;\r
 763           if (clus->GetM20()          > 0.4  ) continue;\r
 764           if (clus->GetM02()          > 0.4  ) continue;\r
 765           if (clus->GetM20()          < 0.03 ) continue;\r
 766           if (clus->GetM02()          < 0.03 ) continue;\r
 767           \r
 768           Double_t x[3];\r
 769           clus->GetPosition(x);\r
 770           TVector3 cluspos(x[0],x[1],x[2]);\r
 771           Double_t deta = teta - cluspos.Eta();\r
 772           Double_t dphi = tphi - cluspos.Phi();\r
 773           if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
 774           if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
 775           fh2dEtadPhi->Fill(deta,dphi);\r
 776           fh2TrackPVsClusterE->Fill(clus->E(),track->P());\r
 777           fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());\r
 778           fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());\r
 779           fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());\r
 780           \r
 781           res = sqrt(dphi*dphi + deta*deta);\r
 782           fh1dR->Fill(res);\r
 783           \r
 784           /////////////////////////////////\r
 785           //Perform electron cut analysis//\r
 786           /////////////////////////////////\r
 787           //Good match\r
 788           if(res < fResidualCut) {\r
 789             fh2dEtadPhiMatched->Fill(deta,dphi);\r
 790             fh2MatchdEdx->Fill(track->P(),dEdx);\r
 791             iCluster = iclus;\r
 792 \r
 793             Double_t energy = clus->E(); \r
 794             if(energy > 0) pOverE = tmom/energy;\r
 795 \r
 796             if (res< minR) {\r
 797               minR  = res;\r
 798               minPe = pOverE;\r
 799               minEp = energy/tmom;\r
 800               minMult = clus->GetNCells() ;\r
 801               minPt = track->Pt();\r
 802             }\r
 803 \r
 804             Int_t cmctag = -1;      \r
 805             if(IsDataMC()) {  \r
 806               //Do you want the cluster or the track label?\r
 807               Int_t input = 0;\r
 808               if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;\r
 809               cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(0),GetReader(),input);\r
 810             }\r
 811             \r
 812             if(fWriteNtuple) {\r
 813               fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,pidProb,imp[0],imp[1]);\r
 814             }\r
 815 \r
 816           } else {\r
 817             fh2dEtadPhiUnmatched->Fill(deta,dphi);\r
 818           }//res cut\r
 819         }//calo cluster loop\r
 820 \r
 821         fh3pOverE->Fill(minPt,minPe ,minMult);\r
 822         fh3EOverp->Fill(minPt,minEp ,minMult);\r
 823         if (trkEle) {\r
 824           fh3pOverE2->Fill(minPt,minPe ,minMult);\r
 825           fh3EOverp2->Fill(minPt,minEp ,minMult);\r
 826         }\r
 827         if (tpcEle) {\r
 828           fh3pOverE3->Fill(minPt,minPe ,minMult);\r
 829           fh3EOverp3->Fill(minPt,minEp ,minMult);\r
 830         }\r
 831         //new\r
 832         if (tmctag>-1 && GetMCSource(tmctag)<8 ) {\r
 833           fh2pOverE->Fill(minPt,minPe );\r
 834           fh2EOverp->Fill(minPt,minEp );\r
 835           if (trkEle) {\r
 836             fh2pOverE2->Fill(minPt,minPe );\r
 837             fh2EOverp2->Fill(minPt,minEp );\r
 838           }\r
 839         }\r
 840 \r
 841         //////////////////////////////\r
 842         //Electron cuts happen here!//\r
 843         //////////////////////////////\r
 844         if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;\r
 845 \r
 846         ///////////////////////////\r
 847         //Fill AOD with electrons//\r
 848         ///////////////////////////\r
 849         //Take all emcal electrons, but the others only if pT < 10 GeV\r
 850         if(emcEle || ( (tpcEle || trkEle) && track->Pt() < 10.) ) {\r
 851 \r
 852           //B-tagging\r
 853           if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");\r
 854           Int_t dvmbtag = GetDVMBtag(track); bt += dvmbtag;\r
 855 \r
 856           fh2EledEdx->Fill(track->P(),dEdx);\r
 857           \r
 858           Double_t eMass = 0.511/1000; //mass in GeV\r
 859           Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);\r
 860           AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);\r
 861           tr.SetLabel(track->GetLabel());\r
 862           tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters\r
 863           tr.SetTrackLabel(itrk,-1); //sets the indices of the original tracks\r
 864           if(emcEle) {//PID determined by EMCAL\r
 865             tr.SetDetector(fCalorimeter);\r
 866           } else {\r
 867             tr.SetDetector("CTS"); //PID determined by CTS\r
 868           }\r
 869 \r
 870           if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);\r
 871           //Make this preserve sign of particle\r
 872           if(track->Charge() < 0) tr.SetPdg(11); //electron is 11\r
 873           else  tr.SetPdg(-11); //positron is -11\r
 874           Int_t btag = 0;\r
 875           if(dvmbtag > 0) tr.SetBTagBit(btag,tr.kDVMTag0);\r
 876           if(dvmbtag > 1) tr.SetBTagBit(btag,tr.kDVMTag1);\r
 877           if(dvmbtag > 2) tr.SetBTagBit(btag,tr.kDVMTag2);\r
 878           tr.SetBtag(btag);\r
 879           \r
 880           //Play with the MC stack if available\r
 881           //Check origin of the candidates\r
 882           if(IsDataMC()){\r
 883             \r
 884             //FIXME:  Need to re-think this for track-oriented analysis\r
 885             //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?\r
 886             tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));\r
 887             \r
 888             if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());\r
 889           }//Work with stack also   \r
 890           \r
 891           AddAODParticle(tr);\r
 892           \r
 893           if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());  \r
 894         }//electron\r
 895       }//pt, fiducial selection\r
 896     }//pid check\r
 897   }//track loop                         \r
 898   \r
 899   fhRefMult->Fill(refmult);\r
 900   fhRefMult2->Fill(refmult2);\r
 901 \r
 902   if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");\r
 903   if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD()  End fill AODs \n");  \r
 904   \r
 905 }\r
 906 \r
 907 //__________________________________________________________________\r
 908 void  AliAnaElectron::MakeAnalysisFillHistograms() \r
 909 {\r
 910   //Do analysis and fill histograms\r
 911 \r
 912   AliStack * stack = 0x0;\r
 913   TParticle * primary = 0x0;\r
 914   AliAODMCParticle * aodprimary = 0x0;\r
 915 \r
 916   Int_t ph1 = 0;  //photonic 1 count\r
 917   Int_t ph2 = 0;  //photonic 2 count\r
 918   Int_t phB = 0;  //both count\r
 919 \r
 920   if(IsDataMC()) {\r
 921     if(GetReader()->ReadStack()){\r
 922       stack =  GetMCStack() ;      \r
 923       if(!stack)\r
 924         printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
 925       \r
 926     }\r
 927   }// is data and MC\r
 928 \r
 929   ////////////////////////////////////\r
 930   //Loop over jets and check for b-tag\r
 931   ////////////////////////////////////\r
 932   Int_t njets = (GetReader()->GetOutputEvent())->GetNJets();\r
 933   if(njets > 0) {\r
 934     if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets.  Performing b-jet tag analysis\n",njets);\r
 935 \r
 936     for(Int_t ijet = 0; ijet < njets ; ijet++) {\r
 937       AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;\r
 938       //Only consider jets with pt > 10 GeV (the rest have to be junk)\r
 939       //printf("AODJet<%d> pt = %2.2f\n",ijet,jet->Pt());\r
 940       if(jet->Pt() < 10.) continue;\r
 941 \r
 942       if(GetDebug() > 3) {\r
 943         printf("AliAODJet ijet = %d\n",ijet);\r
 944         jet->Print("");\r
 945       }\r
 946       //Skip jets not inside a smaller fiducial volume to ensure that\r
 947       //they are completely contained in the EMCAL\r
 948       if(TMath::Abs(jet->Eta()) > fJetEtaCut) continue;\r
 949       if(jet->Phi() < fJetPhiMin || jet->Phi() > fJetPhiMax) continue;\r
 950 \r
 951       //To "tag" the jet, we will look for it to pass our various criteria\r
 952       //For e jet tag, we just look to see which ones have NPEs\r
 953       //For DVM jet tag, we will look for DVM electrons\r
 954       //For IPSig, we compute the IPSig for all tracks and if the\r
 955       //number passing is above the cut, it passes\r
 956       Bool_t leadJet  = kFALSE;\r
 957       if (ijet==0) leadJet= kTRUE;\r
 958       Bool_t eJet = kFALSE;  \r
 959       Bool_t eJet2    = kFALSE; //electron triggered\r
 960       Bool_t hadJet   = kFALSE; //hadron triggered \r
 961       Bool_t dvmJet = kFALSE;  \r
 962       Bool_t ipsigJet = kFALSE;\r
 963       TRefArray* rt = jet->GetRefTracks();\r
 964       Int_t ntrk = rt->GetEntries();\r
 965       Int_t trackCounter[4] = {0,0,0,0}; //for ipsig\r
 966       for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
 967         AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
 968         if( GetIPSignificance(jetTrack, jet->Phi()) > fIPSigCut) trackCounter[0]++;\r
 969         if( GetIPSignificance(jetTrack, jet->Phi()) > 4.) trackCounter[1]++;\r
 970         if( GetIPSignificance(jetTrack, jet->Phi()) > 3.) trackCounter[2]++;\r
 971         if( GetIPSignificance(jetTrack, jet->Phi()) > 2.) trackCounter[3]++;\r
 972         Bool_t isNPE = CheckTrack(jetTrack,"NPE");\r
 973         if(isNPE) eJet = kTRUE;\r
 974         if ( isNPE && jetTrack->Pt()>10.0 ) eJet2  =kTRUE;\r
 975         if (!isNPE && jetTrack->Pt()>10.0) hadJet =kTRUE;\r
 976         Bool_t isDVM = CheckTrack(jetTrack,"DVM");\r
 977         if(isDVM) dvmJet = kTRUE;\r
 978       }\r
 979       fhIPSigBtagQA1->Fill(trackCounter[0]);\r
 980       if(trackCounter[1]>0) fhTagJetPt1x4->Fill(jet->Pt());\r
 981       if(trackCounter[2]>1) fhTagJetPt2x3->Fill(jet->Pt());\r
 982       if(trackCounter[3]>2) fhTagJetPt3x2->Fill(jet->Pt());\r
 983 \r
 984       if(trackCounter[1]>0 && eJet) fhePlusTagJetPt1x4->Fill(jet->Pt());\r
 985       if(trackCounter[2]>1 && eJet) fhePlusTagJetPt2x3->Fill(jet->Pt());\r
 986       if(trackCounter[3]>2 && eJet) fhePlusTagJetPt3x2->Fill(jet->Pt());\r
 987 \r
 988       if(trackCounter[0] > fNTagTrkCut) ipsigJet = kTRUE;\r
 989 \r
 990       if(IsDataMC()) {\r
 991         //determine tagging efficiency & mis-tagging rate\r
 992         //using b-quarks from stack\r
 993         Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());\r
 994         Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());\r
 995         if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");\r
 996         if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");\r
 997         if (dvmJet && GetDebug() > 0)     printf("== found DVM jet==\n");\r
 998 \r
 999         if(isTrueBjet && dvmJet) fhDVMJet->Fill(0.,jet->Pt()); // good tagged\r
1000         if(isTrueBjet && !dvmJet) fhDVMJet->Fill(1.,jet->Pt()); // missed tagged\r
1001         if(!isTrueBjet && dvmJet) fhDVMJet->Fill(2.,jet->Pt());  // fake tagged\r
1002         if(!isTrueBjet && !dvmJet) fhDVMJet->Fill(3.,jet->Pt());  // others\r
1003         if(isTrueDjet && !isTrueBjet &&   dvmJet) fhDVMJet->Fill(4.,jet->Pt()); // charm-tagged\r
1004         if(isTrueDjet && !isTrueBjet &&  !dvmJet) fhDVMJet->Fill(5.,jet->Pt()); // charm -not tagged\r
1005         if(!(isTrueDjet||isTrueBjet ) &&  dvmJet) fhDVMJet->Fill(6.,jet->Pt()); // light flavor -tagged\r
1006         if(!(isTrueDjet||isTrueBjet ) && !dvmJet) fhDVMJet->Fill(7.,jet->Pt()); // light flavor -not tagged\r
1007         if(isTrueBjet &&  eJet && dvmJet) fhDVMJet->Fill(8.,jet->Pt()); // bjet with electron\r
1008         if(isTrueBjet && !eJet && dvmJet) fhDVMJet->Fill(9.,jet->Pt()); // needs more thought\r
1009 \r
1010         if(isTrueBjet) {\r
1011           if(trackCounter[1]>0) fhBJetPt1x4->Fill(jet->Pt());\r
1012           if(trackCounter[2]>1) fhBJetPt2x3->Fill(jet->Pt());\r
1013           if(trackCounter[3]>2) fhBJetPt3x2->Fill(jet->Pt());\r
1014         } else {\r
1015           if(trackCounter[1]>0) fhFakeJetPt1x4->Fill(jet->Pt());\r
1016           if(trackCounter[2]>1) fhFakeJetPt2x3->Fill(jet->Pt());\r
1017           if(trackCounter[3]>2) fhFakeJetPt3x2->Fill(jet->Pt());\r
1018         }\r
1019       }\r
1020 \r
1021       //Fill bjet histograms here\r
1022       if(!(eJet || ipsigJet || dvmJet)) fhJetType->Fill(0.,jet->Pt()); //none\r
1023       if(eJet && !(ipsigJet || dvmJet)) fhJetType->Fill(1.,jet->Pt()); //only ejet\r
1024       if(dvmJet && !(eJet || ipsigJet)) fhJetType->Fill(2.,jet->Pt()); //only dvm\r
1025       if(ipsigJet && !(eJet || dvmJet)) fhJetType->Fill(3.,jet->Pt()); //only ipsig\r
1026       if(eJet && dvmJet && !ipsigJet)   fhJetType->Fill(4.,jet->Pt()); //ejet & dvm\r
1027       if(eJet && ipsigJet && !dvmJet)   fhJetType->Fill(5.,jet->Pt()); //ejet & ipsig\r
1028       if(dvmJet && ipsigJet && !eJet)   fhJetType->Fill(6.,jet->Pt()); //dvm & ipsig\r
1029       if(dvmJet && ipsigJet && eJet)    fhJetType->Fill(7.,jet->Pt()); //all\r
1030       if(dvmJet || ipsigJet || eJet)    fhJetType->Fill(8.,jet->Pt()); //any of them\r
1031       if(eJet  )    fhJetType->Fill(9.,jet->Pt()); //any of them\r
1032       if(dvmJet)    fhJetType->Fill(10.,jet->Pt()); //any of them\r
1033       if(eJet2 )    fhJetType->Fill(11.,jet->Pt()); //any of them\r
1034       if(hadJet)    fhJetType->Fill(12.,jet->Pt()); //any of them\r
1035 \r
1036       if(eJet || ipsigJet || dvmJet) fhBJetEtaPhi->Fill(jet->Eta(),jet->Phi());\r
1037       else fhNonBJetEtaPhi->Fill(jet->Eta(),jet->Phi());\r
1038 \r
1039       //leading jets\r
1040       if (leadJet) {\r
1041         fhLeadJetType->Fill(0.,jet->Pt()); //all\r
1042         if(eJet   )    fhLeadJetType->Fill(1.,jet->Pt());\r
1043         if(eJet2  )    fhLeadJetType->Fill(2.,jet->Pt());\r
1044         if(hadJet )    fhLeadJetType->Fill(3.,jet->Pt());\r
1045         if(eJet   && (dvmJet || ipsigJet) )    fhLeadJetType->Fill(4.,jet->Pt());\r
1046         if(eJet2  && (dvmJet || ipsigJet) )    fhLeadJetType->Fill(5.,jet->Pt());\r
1047         if(hadJet && (dvmJet || ipsigJet) )    fhLeadJetType->Fill(6.,jet->Pt());\r
1048       }\r
1049 \r
1050       for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
1051         AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
1052         Double_t xsi = TMath::Log(jet->Pt()/jetTrack->Pt());\r
1053         if(eJet || ipsigJet || dvmJet) {\r
1054           if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms - We have a bjet!\n");\r
1055           fhBJetXsiFF->Fill(xsi,jet->Pt());\r
1056           fhBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());\r
1057         } else {\r
1058           //Fill non-bjet histograms here\r
1059           fhNonBJetXsiFF->Fill(xsi,jet->Pt());\r
1060           fhNonBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());\r
1061         }\r
1062       }\r
1063 \r
1064     } //jet loop\r
1065   } //jets exist\r
1066   \r
1067   //////////////////////////////\r
1068   //Loop on stored AOD electrons\r
1069   //////////////////////////////\r
1070   Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
1071   if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);\r
1072   \r
1073   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
1074     AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
1075     Int_t pdg = ele->GetPdg();\r
1076     \r
1077     if(GetDebug() > 3) \r
1078       printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;\r
1079     \r
1080     if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue; \r
1081     \r
1082     if(GetDebug() > 1) \r
1083       printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;\r
1084 \r
1085     //MC tag of this electron\r
1086     Int_t mctag = ele->GetTag();\r
1087 \r
1088     //Filter for photonic electrons based on opening angle and Minv\r
1089     //cuts, also fill histograms\r
1090     Bool_t photonic = kFALSE;\r
1091     Bool_t photonic1 = kFALSE;\r
1092     photonic1 = PhotonicPrim(ele); //check against primaries\r
1093     if(photonic1) ph1++;\r
1094     Bool_t photonic2 = kFALSE;\r
1095     photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
1096     if(photonic2) ph2++;\r
1097     if(photonic1 && photonic2) phB++;\r
1098     if(photonic1 || photonic2) photonic = kTRUE;\r
1099 \r
1100     //Fill electron histograms \r
1101     Float_t ptele = ele->Pt();\r
1102     Float_t phiele = ele->Phi();\r
1103     Float_t etaele = ele->Eta();\r
1104 \r
1105     //"Best reconstructed electron spectrum" = EMCAL or tracking\r
1106     //detectors say it is an electron and it does not form a V0\r
1107     //with Minv near a relevant resonance\r
1108     if(!photonic) {\r
1109       fhPtNPEleTTE->Fill(ptele,0); //0 = no MC info\r
1110       if(ele->GetDetector() == fCalorimeter) fhPtNPEleEMCAL->Fill(ptele,0);\r
1111       if(IsDataMC()) {\r
1112         fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag));\r
1113         if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag));\r
1114         if(GetMCSource(mctag) == 1) { //it's a bottom electron, now\r
1115                                       //get the parent's pt\r
1116           Double_t ptbHadron = GetBParentPt(ele->GetLabel());\r
1117           fhPtNPEBHadron->Fill(ptele,ptbHadron);\r
1118         } //mctag\r
1119       } //isdatamc\r
1120     } //!photonic\r
1121 \r
1122     //kept for historical reasons?\r
1123     fhPtElectron  ->Fill(ptele);\r
1124     fhPhiElectron ->Fill(ptele,phiele);\r
1125     fhEtaElectron ->Fill(ptele,etaele);\r
1126 \r
1127     if(photonic) {\r
1128       fhPtPE->Fill(ptele);\r
1129       fhPhiPE->Fill(ptele,phiele);\r
1130       fhEtaPE->Fill(ptele,etaele);\r
1131     } else {\r
1132       fhPtNPE->Fill(ptele);\r
1133       fhPhiNPE->Fill(ptele,phiele);\r
1134       fhEtaNPE->Fill(ptele,etaele);\r
1135     }\r
1136 \r
1137     if(IsDataMC()){\r
1138       if(GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCConversion)){\r
1139         fhPhiConversion ->Fill(ptele,phiele);\r
1140         fhEtaConversion ->Fill(ptele,etaele);\r
1141       }\r
1142     }//Histograms with MC\r
1143     \r
1144   }// aod loop\r
1145 \r
1146   ////////////////////////////////////////////////////////\r
1147   //Fill histograms of pure MC kinematics from the stack//\r
1148   ////////////////////////////////////////////////////////\r
1149   if(IsDataMC()) {\r
1150 \r
1151     //MC Jets\r
1152     TVector3 bjetVect[4];\r
1153     Int_t nPythiaGenJets = 0;\r
1154     AliGenPythiaEventHeader*  pythiaGenHeader = (AliGenPythiaEventHeader*)GetReader()->GetGenEventHeader();\r
1155     if(pythiaGenHeader){\r
1156       //Get Jets from MC header\r
1157       nPythiaGenJets = pythiaGenHeader->NTriggerJets();\r
1158       Int_t iCount = 0;\r
1159       for(int ip = 0;ip < nPythiaGenJets;++ip){\r
1160         if (iCount>3) break;\r
1161         Float_t p[4];\r
1162         pythiaGenHeader->TriggerJet(ip,p);\r
1163         TVector3 tempVect(p[0],p[1],p[2]);\r
1164         if ( TMath::Abs(tempVect.Eta())>fJetEtaCut || tempVect.Phi() < fJetPhiMin || tempVect.Phi() > fJetPhiMax) continue;\r
1165         //Only store it if it has a b-quark within dR < 0.2 of jet axis ?\r
1166         if(IsMcBJet(tempVect.Eta(),tempVect.Phi())) {\r
1167           bjetVect[iCount].SetXYZ(p[0], p[1], p[2]);\r
1168           iCount++;\r
1169         }\r
1170       }\r
1171     }\r
1172 \r
1173     Int_t nPart = GetNumAODMCParticles();\r
1174     if(GetReader()->ReadStack()) nPart = stack->GetNtrack();\r
1175 \r
1176     for(Int_t ipart = 0; ipart < nPart; ipart++) {\r
1177 \r
1178       //All the variables we want from MC particles\r
1179       Double_t px = 0.; Double_t py = 0.; Double_t pz = 0.; Double_t e = 0.;\r
1180       Double_t vx = -999.; Double_t vy = -999.; Double_t vz = -999.; Double_t vt = -999.;\r
1181       Int_t pdg = 0; Int_t mpdg = 0; Double_t mpt = 0.;\r
1182 \r
1183       if(GetReader()->ReadStack()) {\r
1184         primary = stack->Particle(ipart);\r
1185         pdg = primary->GetPdgCode();\r
1186         px = primary->Px(); py = primary->Py(); pz = primary->Pz(); e = primary->Energy();\r
1187         vx = primary->Vx(); vy = primary->Vy(); vz = primary->Vz(); vt = primary->T();\r
1188         if(primary->GetMother(0)>=0) {\r
1189           TParticle *parent = stack->Particle(primary->GetMother(0));\r
1190           if (parent) {\r
1191             mpdg = parent->GetPdgCode();\r
1192             mpt = parent->Pt();\r
1193           }\r
1194         }\r
1195       } else if(GetReader()->ReadAODMCParticles()) {\r
1196         aodprimary =  (AliAODMCParticle*)GetMCParticle(ipart);\r
1197         pdg = aodprimary->GetPdgCode();\r
1198         px = aodprimary->Px(); py = aodprimary->Py(); pz = aodprimary->Pz(); e = aodprimary->E();\r
1199         vx = aodprimary->Xv(); vy = aodprimary->Yv(); vz = aodprimary->Zv(); vt = aodprimary->T();\r
1200         Int_t parentId = aodprimary->GetMother();\r
1201         if(parentId>=0) {\r
1202           AliAODMCParticle *parent = (AliAODMCParticle*)GetMCParticle(parentId);\r
1203           if (parent) {\r
1204             mpdg = parent->GetPdgCode();\r
1205             mpt = parent->Pt();\r
1206           }\r
1207         }       \r
1208       }\r
1209 \r
1210       TLorentzVector mom(px,py,pz,e);\r
1211       TLorentzVector pos(vx,vy,vz,vt);\r
1212       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter);\r
1213       if(mom.Pt() < GetMinPt()) continue;\r
1214       if(!in) continue;\r
1215 \r
1216       if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 2212)\r
1217         fhPtMCHadron->Fill(mom.Pt());\r
1218       \r
1219       //we only care about electrons\r
1220       if(TMath::Abs(pdg) != 11) continue;\r
1221       //we only want TRACKABLE electrons (TPC 85-250cm)\r
1222       if(pos.Rho() > 200.) continue;\r
1223       //Ignore low pt electrons\r
1224       if(mom.Pt() < 0.2) continue;\r
1225         \r
1226       //find out what the ancestry of this electron is\r
1227       Int_t mctag = -1;\r
1228       Int_t input = 0;\r
1229       mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);\r
1230       \r
1231       if(GetMCSource(mctag)==1) { //bottom electron\r
1232         //See if it is within dR < 0.4 of a bjet\r
1233         for(Int_t ij = 0; ij < nPythiaGenJets; ij++) {\r
1234           Double_t deta = primary->Eta() - bjetVect[ij].Eta();\r
1235           Double_t dphi = primary->Phi() - bjetVect[ij].Phi();\r
1236           Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi);\r
1237           if(dR < 0.4) {\r
1238             fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt());\r
1239           }\r
1240         }\r
1241       }\r
1242 \r
1243       if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
1244           (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
1245         {\r
1246           fhMCBHadronElePt->Fill(mom.Pt(), mpt); \r
1247         }\r
1248       //CHECK THAT THIS IS CORRECTLY FILLED - SHOULD WE USE MCSOURCE HERE?\r
1249       fhPtMCElectron->Fill(mom.Pt(),0);  //0 = unfiltered\r
1250       fhPtMCElectron->Fill(mom.Pt(),GetMCSource(mctag));\r
1251 \r
1252       if(GetMCSource(mctag) == 4) {//conversion\r
1253         fhMCXYConversion->Fill(vx,vy);\r
1254         fhMCRadPtConversion->Fill(TMath::Sqrt(vx*vx+vy*vy),mom.Pt());\r
1255       }\r
1256         \r
1257       //fill ntuple\r
1258       if(fWriteNtuple) {\r
1259         fMCEleNtuple->Fill(mctag,mom.Pt(),mom.Phi(),mom.Eta(),vx,vy,vz);\r
1260       }\r
1261     }\r
1262   } //MC loop\r
1263   \r
1264   //if(GetDebug() > 0) \r
1265   printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB);\r
1266 }\r
1267 \r
1268 //__________________________________________________________________\r
1269 Int_t AliAnaElectron::GetDVMBtag(AliAODTrack * tr )\r
1270 {\r
1271   //This method uses the Displaced Vertex between electron-hadron\r
1272   //pairs and the primary vertex to determine whether an electron is\r
1273   //likely from a B hadron.\r
1274 \r
1275   Int_t ncls1 = 0;\r
1276   for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;\r
1277 \r
1278   fhDVMBtagQA3->Fill(ncls1);\r
1279   if (ncls1 < fITSCut) return 0;\r
1280 \r
1281   Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
1282   Bool_t dcaOkay = GetDCA(tr,imp,cov);  //homegrown dca calculation until AOD is fixed                  \r
1283   if(!dcaOkay) {\r
1284     printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d",tr->GetID());\r
1285     return 0;\r
1286   }\r
1287 \r
1288   fhDVMBtagQA4->Fill(imp[0]);\r
1289   if (TMath::Abs(imp[0])   > fImpactCut ) return 0;\r
1290   fhDVMBtagQA5->Fill(imp[1]);\r
1291   if (TMath::Abs(imp[1])   > fImpactCut ) return 0;\r
1292 \r
1293   Int_t nvtx1 = 0;\r
1294   Int_t nvtx2 = 0;\r
1295   Int_t nvtx3 = 0;\r
1296 \r
1297   for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {\r
1298     //loop over assoc\r
1299     AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
1300     Int_t id1 = tr->GetID();\r
1301     Int_t id2 = track2->GetID();\r
1302     if(id1 == id2) continue;\r
1303 \r
1304     Int_t ncls2 = 0;\r
1305     for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;\r
1306     if (ncls2 < fITSCut) continue;\r
1307 \r
1308     if(track2->Pt() < fAssocPtCut) continue;\r
1309 \r
1310     Double_t dphi = tr->Phi() - track2->Phi();\r
1311     if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
1312     if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
1313     Double_t deta = tr->Eta() - track2->Eta();\r
1314     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
1315 \r
1316     if(dr > fDrCut) continue;\r
1317     \r
1318     Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);\r
1319     if (sDca1 > fSdcaCut) nvtx1++;\r
1320     Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);\r
1321     if (sDca2 > fSdcaCut) nvtx2++;\r
1322     Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);\r
1323     if (sDca3 > fSdcaCut) nvtx3++;\r
1324 \r
1325   } //loop over hadrons\r
1326 \r
1327   if(GetDebug() > 0) {\r
1328     if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);\r
1329     if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);\r
1330     if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);\r
1331   }\r
1332 \r
1333   //fill QA histograms\r
1334   fhDVMBtagCut1->Fill(nvtx1,tr->Pt());\r
1335   fhDVMBtagCut2->Fill(nvtx2,tr->Pt());\r
1336   fhDVMBtagCut3->Fill(nvtx3,tr->Pt());\r
1337 \r
1338   return nvtx2;\r
1339 \r
1340 }\r
1341 \r
1342 //__________________________________________________________________\r
1343 Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)\r
1344 {\r
1345   //Compute the signed dca between two tracks\r
1346   //and return the result\r
1347 \r
1348   Double_t signDca=-999.;\r
1349   if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);\r
1350 \r
1351   //=====Now calculate DCA between both tracks=======  \r
1352   Double_t massE = 0.000511;\r
1353   Double_t massK = 0.493677;\r
1354 \r
1355   Double_t bfield = 5.; //kG\r
1356   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
1357 \r
1358   Double_t vertex[3] = {-999.,-999.,-999}; //vertex\r
1359   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
1360     GetReader()->GetVertex(vertex); //If only one file, get the vertex from there\r
1361     //FIXME:  Add a check for whether file 2 is PYTHIA or HIJING\r
1362     //If PYTHIA, then set the vertex from file 2, if not, use the\r
1363     //vertex from file 1\r
1364     if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);\r
1365   }\r
1366   \r
1367   TVector3 primV(vertex[0],vertex[1],vertex[2]) ;\r
1368 \r
1369   if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;\r
1370 \r
1371   AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);\r
1372   AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);\r
1373 \r
1374   Double_t xplane1 = 0.; Double_t xplane2 = 0.;\r
1375   Double_t pairdca = param1->GetDCA(param2,bfield,xplane1,xplane2);\r
1376 \r
1377   param1->PropagateTo(xplane1,bfield);\r
1378   param2->PropagateTo(xplane2,bfield);\r
1379 \r
1380   Int_t id1 = 0, id2 = 0;\r
1381   AliESDv0 bvertex(*param1,id1,*param2,id2);\r
1382   Double_t vx,vy,vz;\r
1383   bvertex.GetXYZ(vx,vy,vz);\r
1384 \r
1385   Double_t emom[3];\r
1386   Double_t hmom[3];\r
1387   param1->PxPyPz(emom);\r
1388   param2->PxPyPz(hmom);\r
1389   TVector3 emomAtB(emom[0],emom[1],emom[2]);\r
1390   TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);\r
1391   TVector3 secvtxpt(vx,vy,vz);\r
1392   TVector3 decayvector(0,0,0);\r
1393   decayvector = secvtxpt - primV; //decay vector from PrimVtx\r
1394   Double_t decaylength = decayvector.Mag();\r
1395 \r
1396   if(GetDebug() > 0) {\r
1397     printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );\r
1398     printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );\r
1399   }\r
1400 \r
1401   if (masscut<1.1) fhDVMBtagQA1->Fill(pairdca,decaylength);\r
1402 \r
1403   if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {\r
1404     TVector3 sumMom = emomAtB+hmomAtB;\r
1405     Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);\r
1406     Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);\r
1407     Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);\r
1408     Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
1409     Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));\r
1410     Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();\r
1411 \r
1412     if (masscut<1.1) fhDVMBtagQA2->Fill(sDca, mass);\r
1413 \r
1414     if (mass > masscut && massPhot > 0.1) signDca = sDca;\r
1415     \r
1416     if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);\r
1417     if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);\r
1418   }\r
1419 \r
1420   //clean up\r
1421   delete param1;\r
1422   delete param2;\r
1423 \r
1424   return signDca;\r
1425 }\r
1426 \r
1427 //__________________________________________________________________\r
1428 Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi)\r
1429 {\r
1430   //get signed impact parameter significance of the given AOD track\r
1431   //for the given jet\r
1432 \r
1433   Int_t trackIndex = 0;\r
1434   Int_t ntrk = GetAODCTS()->GetEntriesFast();\r
1435   for (Int_t k2 =0; k2 < ntrk ; k2++) {\r
1436     //loop over assoc\r
1437     AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
1438     int id1 = tr->GetID();\r
1439     int id2 = track2->GetID();\r
1440     if(id1 == id2) {\r
1441       trackIndex = k2;//FIXME: check if GetAODCTS stores tracks in the\r
1442                       //same order of the event\r
1443       break;\r
1444     }\r
1445   }\r
1446 \r
1447   Double_t significance=0;\r
1448   Double_t magField = 0;\r
1449   Double_t maxD = 10000.;\r
1450   Double_t impPar[] = {0,0};\r
1451   Double_t ipCov[]={0,0,0};\r
1452   Double_t ipVec2D[] = {0,0};\r
1453 \r
1454   AliVEvent* vEvent = (AliVEvent*)GetReader()->GetInputEvent();\r
1455   if(!vEvent) return -97;\r
1456   AliVVertex* vv = (AliVVertex*)vEvent->GetPrimaryVertex();\r
1457   if(!vv) return -98;\r
1458   AliVTrack* vTrack = (AliVTrack*)vEvent->GetTrack(trackIndex);\r
1459   if(!vTrack) return -99;\r
1460   AliESDtrack esdTrack(vTrack);\r
1461   if(!esdTrack.PropagateToDCA(vv, magField, maxD, impPar, ipCov)) return -100;\r
1462   if(ipCov[0]<0) return -101;\r
1463 \r
1464   Double_t Pxy[] = {esdTrack.Px(), esdTrack.Py()};\r
1465   Double_t Txy[] = {esdTrack.Xv(), esdTrack.Yv()};\r
1466   Double_t Vxy[] = {vv->GetX(),  vv->GetY()};\r
1467   GetImpactParamVect(Pxy, Txy, Vxy, ipVec2D);\r
1468         Double_t phiIP = TMath::ATan2(ipVec2D[1], ipVec2D[0]) + (TMath::Abs(ipVec2D[1])-ipVec2D[1])/TMath::Abs(ipVec2D[1])*TMath::Pi();\r
1469   Double_t cosTheta = TMath::Cos(jetPhi - phiIP);\r
1470   Double_t sign = cosTheta/TMath::Abs(cosTheta);\r
1471   significance = TMath::Abs(impPar[0])/TMath::Sqrt(ipCov[0])*sign;\r
1472   //ip = fabs(impPar[0]);\r
1473   fhIPSigBtagQA2->Fill(significance);\r
1474   return significance;\r
1475 }\r
1476 \r
1477 //__________________________________________________________________\r
1478 void AliAnaElectron::GetImpactParamVect(Double_t Pxy[2], Double_t Txy[2], Double_t Vxy[2], Double_t IPxy[2])\r
1479 {\r
1480   //px,py: momentum components at the origin of the track; tx, ty:\r
1481   //origin (x,y) of track; vx, vy: coordinates of primary vertex\r
1482   // analytical geometry auxiliary variables\r
1483   Double_t mr = Pxy[1]/Pxy[0]; //angular coeficient of the straight\r
1484                               //line that lies on top of track\r
1485                               //momentum\r
1486   Double_t br = Txy[1] - mr*Txy[0]; //linear coeficient of the straight\r
1487                                    //line that lies on top of track\r
1488                                    //momentum\r
1489   Double_t ms = -1./mr; //angular coeficient of the straight line that\r
1490                        //lies on top of the impact parameter line\r
1491   //  Double_t bs = Vxy[1] - ms*Vxy[0]; //linear coeficient of the straight\r
1492                                    //line that lies on top of the\r
1493                                    //impact parameter line \r
1494   Double_t xIntersection = (mr*Txy[0] - ms*Vxy[0] + Vxy[1] - Txy[1])/(mr - ms);\r
1495   Double_t yIntersection = mr*xIntersection + br;\r
1496   //if(ceil(10000*yIntersection) - ceil(10000*(ms*xIntersection + bs))\r
1497   //!= 0 )cout<<yIntersection<<", "<<ms*xIntersection + bs<<endl;\r
1498   IPxy[0] = xIntersection - Vxy[0];\r
1499   IPxy[1] = yIntersection - Vxy[1];\r
1500   return;\r
1501 }\r
1502 \r
1503 //__________________________________________________________________\r
1504 Bool_t AliAnaElectron::PhotonicPrim(const AliAODPWG4Particle* part) \r
1505 {\r
1506   //This method checks the opening angle and invariant mass of\r
1507   //electron pairs within the AliAODPWG4Particle list to see if \r
1508   //they are likely to be photonic electrons\r
1509 \r
1510   Bool_t itIS = kFALSE;\r
1511 \r
1512   Double_t massE = 0.000511;\r
1513   Double_t massEta = 0.547;\r
1514   Double_t massRho0 = 0.770;\r
1515   Double_t massOmega = 0.782;\r
1516   Double_t massPhi = 1.020;\r
1517 \r
1518   Double_t bfield = 5.; //kG\r
1519   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
1520 \r
1521   Int_t pdg1 = part->GetPdg();\r
1522   Int_t trackId = part->GetTrackLabel(0);\r
1523   AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);\r
1524   if(!track) {\r
1525     if(GetDebug() > 0) printf("AliAnaElectron::PhotonicPrim - can't get the AOD Track from the particle!  Skipping the photonic check");\r
1526     return kFALSE; //Don't proceed because we can't get the track\r
1527   }\r
1528 \r
1529   AliExternalTrackParam *param1 = new AliExternalTrackParam(track);\r
1530 \r
1531   //Loop on stored AOD electrons and compute the angle differences and Minv\r
1532   for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {\r
1533     AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);\r
1534     Int_t track2Id = part2->GetTrackLabel(0);\r
1535     if(trackId == track2Id) continue;\r
1536     Int_t pdg2 = part2->GetPdg();\r
1537     if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;\r
1538     if(part2->GetDetector() != fCalorimeter) continue;\r
1539 \r
1540     //JLK: Check opp. sign pairs only\r
1541     if(pdg1*pdg2 > 0) continue; //skip same-sign pairs\r
1542 \r
1543     //propagate to common vertex and check opening angle\r
1544     AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(track2Id);\r
1545     if(!track2) {\r
1546       if(GetDebug() >0) printf("AliAnaElectron::PhotonicPrim - problem getting the partner track.  Continuing on to the next one");\r
1547       continue;\r
1548     }\r
1549     AliExternalTrackParam *param2 = new AliExternalTrackParam(track2);\r
1550     Int_t id1 = 0, id2 = 0;\r
1551     AliESDv0 photonVtx(*param1,id1,*param2,id2);\r
1552     Double_t vx,vy,vz;\r
1553     photonVtx.GetXYZ(vx,vy,vz);\r
1554 \r
1555     Double_t p1mom[3];\r
1556     Double_t p2mom[3];\r
1557     param1->PxPyPz(p1mom);\r
1558     param2->PxPyPz(p2mom);\r
1559 \r
1560     TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);\r
1561     TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);\r
1562     TVector3 sumMom = p1momAtB+p2momAtB;\r
1563 \r
1564     Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);\r
1565     Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);\r
1566     Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
1567 \r
1568     Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);\r
1569     fh1OpeningAngle->Fill(dphi);\r
1570     fh1MinvPhoton->Fill(mass);\r
1571 \r
1572     if(mass < 0.1 ||\r
1573        (mass > massEta-0.05 || mass < massEta+0.05) ||\r
1574        (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
1575        (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
1576        (mass > massPhi-0.05 || mass < massPhi+0.05)) \r
1577       {\r
1578       \r
1579         if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");\r
1580         itIS = kTRUE;\r
1581       }\r
1582     \r
1583     //clean up\r
1584     delete param2;\r
1585     \r
1586   }\r
1587 \r
1588   delete param1;\r
1589   return itIS;\r
1590 \r
1591 }\r
1592 \r
1593 //__________________________________________________________________\r
1594 Bool_t AliAnaElectron::PhotonicV0(Int_t id) \r
1595 {\r
1596   //This method checks to see whether a track that has been flagged as\r
1597   //an electron was determined to match to a V0 candidate with\r
1598   //invariant mass consistent with photon conversion\r
1599 \r
1600   Bool_t itIS = kFALSE;\r
1601 \r
1602   Double_t massEta = 0.547;\r
1603   Double_t massRho0 = 0.770;\r
1604   Double_t massOmega = 0.782;\r
1605   Double_t massPhi = 1.020;\r
1606   \r
1607   //---Get V0s---\r
1608   AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();\r
1609   int nv0s = aod->GetNumberOfV0s();\r
1610   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
1611     AliAODv0 *v0 = aod->GetV0(iV0);\r
1612     if (!v0) continue;\r
1613     double radius = v0->RadiusV0();\r
1614     double mass = v0->InvMass2Prongs(0,1,11,11);\r
1615     if(GetDebug() > 0) {\r
1616       printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );\r
1617       printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);\r
1618       printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );\r
1619     }\r
1620     if (mass < 0.100 ||\r
1621         (mass > massEta-0.05 || mass < massEta+0.05) ||\r
1622         (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
1623         (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
1624         (mass > massPhi-0.05 || mass < massPhi+0.05)) {\r
1625       if ( id == v0->GetNegID() || id == v0->GetPosID()) {\r
1626         itIS=kTRUE;\r
1627         if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );\r
1628       }\r
1629     } }\r
1630   return itIS;\r
1631 \r
1632 }\r
1633 \r
1634 //__________________________________________________________________\r
1635 Bool_t AliAnaElectron::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3]) \r
1636 {\r
1637   //Use the Event vertex and AOD track information to get\r
1638   //a real impact parameter for the track\r
1639   //Once alice-off gets its act together and fixes the AOD, this\r
1640   //should become obsolete.\r
1641 \r
1642   Double_t bfield = 5.; //kG\r
1643   Double_t maxD = 100000.; //max transverse IP\r
1644   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
1645     bfield = GetReader()->GetBField();\r
1646     AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();\r
1647     AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();\r
1648     AliESDtrack esdTrack(track);\r
1649     Bool_t gotit = esdTrack.PropagateToDCA(vv,bfield,maxD,impPar,cov);\r
1650     return gotit;\r
1651   }\r
1652 \r
1653   return kFALSE;\r
1654 \r
1655 }\r
1656 \r
1657 //__________________________________________________________________\r
1658 Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type) \r
1659 {\r
1660   //Check this track to see if it is also tagged as an electron in the\r
1661   //AliAODPWG4Particle list and if it is non-photonic\r
1662 \r
1663   Bool_t pass = kFALSE;\r
1664 \r
1665   Int_t trackId = track->GetID(); //get the index in the reader\r
1666 \r
1667   Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
1668   if(GetDebug() > 3) printf("AliAnaElectron::CheckTrack() - aod branch entries %d\n", naod);\r
1669   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
1670     AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
1671     Int_t label = ele->GetTrackLabel(0);\r
1672     if(label != trackId) continue;  //skip to the next one if they don't match\r
1673 \r
1674     if(strcmp(type,"DVM")==0) { \r
1675       if(ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag1) ||\r
1676          ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag2))\r
1677         pass = kTRUE;\r
1678 \r
1679     } else if (strcmp(type,"NPE")==0) {\r
1680 \r
1681       Bool_t photonic = kFALSE;\r
1682       Bool_t photonic1 = kFALSE;\r
1683       photonic1 = PhotonicPrim(ele); //check against primaries\r
1684       Bool_t photonic2 = kFALSE;\r
1685       photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
1686       if(photonic1 || photonic2) photonic = kTRUE;\r
1687       \r
1688       if(!photonic) pass = kTRUE;\r
1689 \r
1690     } else {\r
1691       return kFALSE;\r
1692     }\r
1693   }\r
1694 \r
1695   return pass;\r
1696 \r
1697 }\r
1698 \r
1699 //__________________________________________________________________\r
1700 Double_t AliAnaElectron::GetBParentPt(Int_t ipart)\r
1701 {\r
1702   //return MC B parent pt\r
1703   if(GetReader()->ReadStack()) { //only done if we have the stack                                                                                               \r
1704     AliStack* stack = GetMCStack();\r
1705     if(!stack) {\r
1706       printf("Problem getting stack\n");\r
1707       return 0.;\r
1708     }\r
1709     TParticle* prim = stack->Particle(ipart);\r
1710     if(prim->GetMother(0)>=0) {\r
1711       Int_t mpdg = 0;\r
1712       TParticle *parent = stack->Particle(prim->GetMother(0));\r
1713       if(parent) mpdg = parent->GetPdgCode();\r
1714 \r
1715       if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
1716           (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
1717         return parent->Pt();\r
1718     }\r
1719   } else if(GetReader()->ReadAODMCParticles()){\r
1720     AliAODMCParticle* prim = (AliAODMCParticle*)GetMCParticle(ipart);\r
1721     if(prim->GetMother()>=0) {\r
1722       Int_t mpdg = 0;\r
1723       AliAODMCParticle* parent = (AliAODMCParticle*)GetMCParticle(prim->GetMother());\r
1724       if(parent) mpdg = parent->GetPdgCode();\r
1725       if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
1726           (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
1727         return parent->Pt();\r
1728     }\r
1729   }\r
1730   return 0.;\r
1731 }\r
1732 \r
1733 //__________________________________________________________________\r
1734 Int_t AliAnaElectron::GetMCSource(Int_t tag)\r
1735 {\r
1736   //For determining how to classify electrons using MC info\r
1737   //the number returned is the bin along one axis of 2-d histograms in\r
1738   //which to fill this electron\r
1739 \r
1740   //Do this first\r
1741   if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;\r
1742 \r
1743   if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) {\r
1744     //Bottom\r
1745     if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 1;\r
1746     //Charm only\r
1747     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC)\r
1748             && !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 2;\r
1749     //Charm from bottom\r
1750     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB)) return 3;\r
1751     //    //Conversion\r
1752     //else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;\r
1753     //Dalitz\r
1754     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) \r
1755        || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) \r
1756        || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) return 5; \r
1757     //W,Z\r
1758     else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay)\r
1759             || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay)) return 6;\r
1760     //Everything else\r
1761     else \r
1762       return 7;\r
1763   } else {\r
1764     //Misidentified electron\r
1765     return 8;\r
1766   }\r
1767 \r
1768 }\r
1769 \r
1770 //__________________________________________________________________\r
1771 Int_t AliAnaElectron::GetNumAODMCParticles() \r
1772 {\r
1773   //Get the number of AliAODMCParticles, if any\r
1774   Int_t num = 0;\r
1775 \r
1776   TClonesArray * mcparticles0 = 0x0;\r
1777   TClonesArray * mcparticles1 = 0x0;\r
1778 \r
1779   if(GetReader()->ReadAODMCParticles()){\r
1780     //Get the list of MC particles\r
1781     //                                                                                                 \r
1782     mcparticles0 = GetReader()->GetAODMCParticles(0);\r
1783     if(!mcparticles0 && GetDebug() > 0) {\r
1784       printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
1785     }\r
1786     if(GetReader()->GetSecondInputAODTree()){\r
1787       mcparticles1 = GetReader()->GetAODMCParticles(1);\r
1788       if(!mcparticles1 && GetDebug() > 0) {\r
1789         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");\r
1790       }\r
1791     }\r
1792 \r
1793     Int_t npart0 = mcparticles0->GetEntriesFast();\r
1794     Int_t npart1 = 0;\r
1795     if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
1796     Int_t npart = npart0+npart1;\r
1797     return npart;\r
1798 \r
1799   }\r
1800 \r
1801   return num;\r
1802 }\r
1803 //__________________________________________________________________\r
1804 AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t ipart) \r
1805 {\r
1806   //Get the MC particle at position ipart\r
1807 \r
1808   AliAODMCParticle* aodprimary = 0x0;\r
1809   TClonesArray * mcparticles0 = 0x0;\r
1810   TClonesArray * mcparticles1 = 0x0;\r
1811 \r
1812   if(GetReader()->ReadAODMCParticles()){\r
1813     //Get the list of MC particles                                                                                                                           \r
1814     mcparticles0 = GetReader()->GetAODMCParticles(0);\r
1815     if(!mcparticles0 && GetDebug() > 0) {\r
1816       printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
1817     }\r
1818     if(GetReader()->GetSecondInputAODTree()){\r
1819       mcparticles1 = GetReader()->GetAODMCParticles(1);\r
1820       if(!mcparticles1 && GetDebug() > 0) {\r
1821         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");\r
1822       }\r
1823     }\r
1824 \r
1825     Int_t npart0 = mcparticles0->GetEntriesFast();\r
1826     Int_t npart1 = 0;\r
1827     if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
1828     if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
1829     else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);\r
1830     if(!aodprimary) {\r
1831       printf("AliAnaElectron::GetMCParticle() *** no primary ***:  label %d \n", ipart);\r
1832       return 0x0;\r
1833     }\r
1834 \r
1835   } else {\r
1836     printf("AliAnaElectron::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
1837   }\r
1838   return aodprimary;\r
1839 \r
1840 }\r
1841 \r
1842 //__________________________________________________________________\r
1843 Bool_t  AliAnaElectron::IsMcBJet(Double_t jeta, Double_t jphi)\r
1844 {\r
1845   //Check the jet eta,phi against that of the b-quark\r
1846   //to decide whether it is an MC B-jet\r
1847   Bool_t bjet=kFALSE;\r
1848 \r
1849   //      printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() );\r
1850 \r
1851   AliStack* stack = 0x0;\r
1852   \r
1853   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
1854 \r
1855     Double_t pphi = -999.;\r
1856     Double_t peta = -999.;\r
1857     Int_t pdg = 0;\r
1858     if(GetReader()->ReadStack()) {\r
1859       stack = GetMCStack();\r
1860       if(!stack) {\r
1861         printf("AliAnaElectron::IsMCBJet() *** no stack ***: \n");\r
1862         return kFALSE;\r
1863       }\r
1864       TParticle* primary = stack->Particle(ipart);\r
1865       if (!primary) continue;\r
1866       pdg = primary->GetPdgCode();\r
1867       pphi = primary->Phi();\r
1868       peta = primary->Eta();\r
1869     } else if(GetReader()->ReadAODMCParticles()) {\r
1870       AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
1871       if(!aodprimary) continue;\r
1872       pdg = aodprimary->GetPdgCode();\r
1873       pphi = aodprimary->Phi();\r
1874       peta = aodprimary->Eta();\r
1875     }\r
1876     if ( TMath::Abs(pdg) != 5) continue;\r
1877     \r
1878     //      printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());\r
1879     Double_t dphi = jphi - pphi;\r
1880     Double_t deta = jeta - peta;\r
1881     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
1882     \r
1883     if (dr < 0.2) {\r
1884       bjet=kTRUE;\r
1885       //printf("MTH: **** found matching MC-Bjet: PDG=%d, pt=%f,dr=%f \n", pdgcode, primary->Pt(),dr );\r
1886       break;\r
1887     }\r
1888   }\r
1889   return bjet;\r
1890 \r
1891 }\r
1892 \r
1893 //__________________________________________________________________\r
1894 Bool_t  AliAnaElectron::IsMcDJet(Double_t jeta, Double_t jphi)\r
1895 {\r
1896   //Check if this jet is a charm jet\r
1897   Bool_t cjet=kFALSE;\r
1898 \r
1899   AliStack* stack = 0x0;\r
1900 \r
1901   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
1902     \r
1903     Double_t pphi = -999.;\r
1904     Double_t peta = -999.;\r
1905     Int_t pdg = 0;\r
1906     if(GetReader()->ReadStack()) {\r
1907       stack = GetMCStack();\r
1908       if(!stack) {\r
1909         printf("AliAnaElectron::IsMCDJet() *** no stack ***: \n");\r
1910         return kFALSE;\r
1911       }\r
1912       TParticle* primary = stack->Particle(ipart);\r
1913       if (!primary) continue;\r
1914       pdg = primary->GetPdgCode();\r
1915       pphi = primary->Phi();\r
1916       peta = primary->Eta();\r
1917     } else if(GetReader()->ReadAODMCParticles()) {\r
1918       AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
1919       if(!aodprimary) continue;\r
1920       pdg = aodprimary->GetPdgCode();\r
1921       pphi = aodprimary->Phi();\r
1922       peta = aodprimary->Eta();\r
1923     }\r
1924 \r
1925     if ( TMath::Abs(pdg) != 4) continue;\r
1926 \r
1927     Double_t dphi = jphi - pphi;\r
1928     Double_t deta = jeta - peta;\r
1929     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
1930     \r
1931     if (dr < 0.2) {\r
1932       cjet=kTRUE;\r
1933       break;\r
1934     }\r
1935   }\r
1936 \r
1937   return cjet;\r
1938 \r
1939 }\r
1940 \r
1941 //__________________________________________________________________\r
1942 void AliAnaElectron::Print(const Option_t * opt) const\r
1943 {\r
1944   //Print some relevant parameters set for the analysis\r
1945   \r
1946   if(! opt)\r
1947     return;\r
1948   \r
1949   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
1950   AliAnaPartCorrBaseClass::Print(" ");\r
1951 \r
1952   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;\r
1953   printf("pOverE range           =     %f - %f\n",fpOverEmin,fpOverEmax);\r
1954   printf("residual cut           =     %f\n",fResidualCut);\r
1955   printf("---DVM Btagging\n");\r
1956   printf("max IP-cut (e,h)       =     %f\n",fImpactCut);\r
1957   printf("min ITS-hits           =     %d\n",fITSCut);\r
1958   printf("max dR (e,h)           =     %f\n",fDrCut);\r
1959   printf("max pairDCA            =     %f\n",fPairDcaCut);\r
1960   printf("max decaylength        =     %f\n",fDecayLenCut);\r
1961   printf("min Associated Pt      =     %f\n",fAssocPtCut);\r
1962   printf("---IPSig Btagging\n");\r
1963   printf("min tag track          =     %d\n",fNTagTrkCut);\r
1964   printf("min IP significance    =     %f\n",fIPSigCut);\r
1965   printf("    \n") ;\r
1966         \r
1967 } \r
1968 \r
1969 //________________________________________________________________________\r
1970 void AliAnaElectron::ReadHistograms(TList* /* outputList */)\r
1971 {\r
1972   // Needed when Terminate is executed in distributed environment                             \r
1973   // Refill analysis histograms of this class with corresponding\r
1974   // histograms in output list.   \r
1975 \r
1976   // Histograms of this analsys are kept in the same list as other\r
1977   // analysis, recover the position of\r
1978   // the first one and then add the next                                                      \r
1979   //Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));\r
1980 \r
1981   //Read histograms, must be in the same order as in\r
1982   //GetCreateOutputObject.                   \r
1983   //fh1pOverE     = (TH1F *) outputList->At(index);\r
1984   //fh1dR         = (TH1F *) outputList->At(index++);\r
1985   //fh2EledEdx    = (TH2F *) outputList->At(index++);\r
1986   //fh2MatchdEdx  = (TH2F *) outputList->At(index++);\r
1987   \r
1988 }\r
1989 \r
1990 //__________________________________________________________________\r
1991 void  AliAnaElectron::Terminate(TList* outputList)\r
1992 {\r
1993 \r
1994   //Do some plots to end\r
1995   //Recover histograms from output histograms list, needed for\r
1996   //distributed analysis.                \r
1997   //ReadHistograms(outputList);\r
1998 \r
1999   printf(" AliAnaElectron::Terminate()  *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;\r
2000 \r
2001 }\r
2002 \r