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