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