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