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