1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\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
17 //_________________________________________________________________________
\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
23 // -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)
\r
24 //////////////////////////////////////////////////////////////////////////////
\r
26 // --- ROOT system ---
\r
29 #include <TParticle.h>
\r
30 #include <TNtuple.h>
\r
31 #include <TClonesArray.h>
\r
32 #include <TObjString.h>
\r
33 //#include <Riostream.h>
\r
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 "AliFiducialCut.h"
\r
41 #include "AliAODTrack.h"
\r
42 #include "AliAODPid.h"
\r
43 #include "AliCaloPID.h"
\r
44 #include "AliAODMCParticle.h"
\r
45 #include "AliStack.h"
\r
46 #include "AliExternalTrackParam.h"
\r
47 #include "AliESDv0.h"
\r
48 #include "AliESDtrack.h"
\r
49 #include "AliAODJet.h"
\r
50 #include "AliAODEvent.h"
\r
51 #include "AliGenPythiaEventHeader.h"
\r
53 ClassImp(AliAnaElectron)
\r
55 //____________________________________________________________________________
\r
56 AliAnaElectron::AliAnaElectron()
\r
57 : AliAnaPartCorrBaseClass(),fCalorimeter(""),
\r
58 fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.),
\r
59 fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),
\r
60 fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),
\r
61 fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),
\r
62 fWriteNtuple(kFALSE),
\r
64 fhImpactXY(0),fhRefMult(0),fhRefMult2(0),
\r
66 fh3pOverE(0),fh3EOverp(0),fh3pOverE2(0),fh3EOverp2(0),fh3pOverE3(0),fh3EOverp3(0),
\r
67 fh2pOverE(0),fh2EOverp(0),fh2pOverE2(0),fh2EOverp2(0),
\r
68 fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),
\r
69 fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),fh2TrackPVsClusterE(0),
\r
70 fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),
\r
71 //Photonic electron checks
\r
72 fh1OpeningAngle(0),fh1MinvPhoton(0),
\r
73 //Reconstructed electrons
\r
74 fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
\r
75 fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),
\r
76 fhPtPE(0),fhPhiPE(0),fhEtaPE(0),
\r
77 //for comparisons with tracking detectors
\r
78 fhPtHadron(0),fhPtNPEleTPC(0),fhPtNPEleTPCTRD(0),fhPtNPEleTTE(0),
\r
81 fhDVMBtagCut1(0),fhDVMBtagCut2(0),fhDVMBtagCut3(0),fhDVMBtagQA1(0),fhDVMBtagQA2(0),
\r
82 fhDVMBtagQA3(0),fhDVMBtagQA4(0),fhDVMBtagQA5(0),
\r
84 fhIPSigBtagQA1(0),fhIPSigBtagQA2(0),fhTagJetPt1x4(0),fhTagJetPt2x3(0),fhTagJetPt3x2(0),
\r
85 fhePlusTagJetPt1x4(0),fhePlusTagJetPt2x3(0),fhePlusTagJetPt3x2(0),
\r
87 fhJetType(0),fhLeadJetType(0),fhBJetXsiFF(0),fhBJetPtFF(0),fhBJetEtaPhi(0),
\r
88 fhNonBJetXsiFF(0),fhNonBJetPtFF(0),fhNonBJetEtaPhi(0),
\r
89 /////////////////////////////////////////////////////////////
\r
90 //Histograms that rely on MC info (not filled for real data)
\r
92 //reco electrons from various sources
\r
93 fhPhiConversion(0),fhEtaConversion(0),
\r
94 //for comparisons with tracking detectors
\r
97 //for computing efficiency of B-jet tags
\r
98 fhBJetPt1x4(0),fhBJetPt2x3(0),fhBJetPt3x2(0),
\r
99 fhFakeJetPt1x4(0),fhFakeJetPt2x3(0),fhFakeJetPt3x2(0),fhDVMJet(0),
\r
100 //MC rate histograms/ntuple
\r
101 fMCEleNtuple(0),fhMCBJetElePt(0),fhMCBHadronElePt(0),fhPtMCHadron(0),fhPtMCElectron(0),
\r
102 fhMCXYConversion(0),fhMCRadPtConversion(0)
\r
106 //Initialize parameters
\r
111 //____________________________________________________________________________
\r
112 AliAnaElectron::AliAnaElectron(const AliAnaElectron & g)
\r
113 : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),
\r
114 fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),
\r
115 fResidualCut(g.fResidualCut),fMinClusEne(g.fMinClusEne),
\r
116 fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),
\r
117 fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),
\r
118 fNTagTrkCut(g.fNTagTrkCut),fIPSigCut(g.fIPSigCut),
\r
119 fJetEtaCut(g.fJetEtaCut),fJetPhiMin(g.fJetPhiMin),fJetPhiMax(g.fJetPhiMax),
\r
120 fWriteNtuple(g.fWriteNtuple),
\r
122 fhImpactXY(g.fhImpactXY),fhRefMult(g.fhRefMult),fhRefMult2(g.fhRefMult2),
\r
124 fh3pOverE(g.fh3pOverE),fh3EOverp(g.fh3EOverp),
\r
125 fh3pOverE2(g.fh3pOverE2),fh3EOverp2(g.fh3EOverp2),
\r
126 fh3pOverE3(g.fh3pOverE3),fh3EOverp3(g.fh3EOverp3),
\r
127 fh2pOverE(g.fh2pOverE),fh2EOverp(g.fh2EOverp),
\r
128 fh2pOverE2(g.fh2pOverE2),fh2EOverp2(g.fh2EOverp2),
\r
129 fh1dR(g.fh1dR),fh2EledEdx(g.fh2EledEdx),
\r
130 fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),
\r
131 fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),
\r
132 fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),
\r
133 fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),
\r
134 //Photonic electron checks
\r
135 fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),
\r
136 //Reconstructed electrons
\r
137 fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron),
\r
138 fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE),
\r
139 fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE),
\r
140 //for comparisons with tracking detectors
\r
141 fhPtHadron(g.fhPtHadron),fhPtNPEleTPC(g.fhPtNPEleTPC),
\r
142 fhPtNPEleTPCTRD(g.fhPtNPEleTPCTRD),fhPtNPEleTTE(g.fhPtNPEleTTE),
\r
143 fhPtNPEleEMCAL(g.fhPtNPEleEMCAL),
\r
145 fhDVMBtagCut1(g.fhDVMBtagCut1),fhDVMBtagCut2(g.fhDVMBtagCut2),fhDVMBtagCut3(g.fhDVMBtagCut3),
\r
146 fhDVMBtagQA1(g.fhDVMBtagQA1),fhDVMBtagQA2(g.fhDVMBtagQA2),
\r
147 fhDVMBtagQA3(g.fhDVMBtagQA3),fhDVMBtagQA4(g.fhDVMBtagQA4),fhDVMBtagQA5(g.fhDVMBtagQA5),
\r
149 fhIPSigBtagQA1(g.fhIPSigBtagQA1),fhIPSigBtagQA2(g.fhIPSigBtagQA2),
\r
150 fhTagJetPt1x4(g.fhTagJetPt1x4),fhTagJetPt2x3(g.fhTagJetPt2x3),fhTagJetPt3x2(g.fhTagJetPt3x2),
\r
151 fhePlusTagJetPt1x4(g.fhePlusTagJetPt1x4),fhePlusTagJetPt2x3(g.fhePlusTagJetPt2x3),
\r
152 fhePlusTagJetPt3x2(g.fhePlusTagJetPt3x2),
\r
154 fhJetType(g.fhJetType),fhLeadJetType(g.fhLeadJetType),fhBJetXsiFF(g.fhBJetXsiFF),
\r
155 fhBJetPtFF(g.fhBJetPtFF),fhBJetEtaPhi(g.fhBJetEtaPhi),fhNonBJetXsiFF(g.fhNonBJetXsiFF),
\r
156 fhNonBJetPtFF(g.fhNonBJetPtFF),fhNonBJetEtaPhi(g.fhNonBJetEtaPhi),
\r
157 /////////////////////////////////////////////////////////////
\r
158 //Histograms that rely on MC info (not filled for real data)
\r
159 fEleNtuple(g.fEleNtuple),
\r
160 //reco electrons from various sources
\r
161 fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
\r
162 //for comparisons with tracking detectors
\r
163 fhPtTrack(g.fhPtTrack),
\r
164 fhPtNPEBHadron(g.fhPtNPEBHadron),
\r
165 //for computing efficiency of B-jet tags
\r
166 fhBJetPt1x4(g.fhBJetPt1x4),fhBJetPt2x3(g.fhBJetPt2x3),
\r
167 fhBJetPt3x2(g.fhBJetPt3x2),
\r
168 fhFakeJetPt1x4(g.fhFakeJetPt1x4),fhFakeJetPt2x3(g.fhBJetPt2x3),
\r
169 fhFakeJetPt3x2(g.fhFakeJetPt3x2),fhDVMJet(g.fhDVMJet),
\r
170 //MC rate histograms/ntuple
\r
171 fMCEleNtuple(g.fMCEleNtuple),fhMCBJetElePt(g.fhMCBJetElePt),
\r
172 fhMCBHadronElePt(g.fhMCBHadronElePt),
\r
173 fhPtMCHadron(g.fhPtMCHadron),fhPtMCElectron(g.fhPtMCElectron),
\r
174 fhMCXYConversion(g.fhMCXYConversion),fhMCRadPtConversion(g.fhMCRadPtConversion)
\r
180 //_________________________________________________________________________
\r
181 AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)
\r
183 // assignment operator
\r
185 if(&g == this) return *this;
\r
186 fCalorimeter = g.fCalorimeter;
\r
187 fpOverEmin = g.fpOverEmin;
\r
188 fpOverEmax = g.fpOverEmax;
\r
189 fResidualCut = g.fResidualCut;
\r
190 fMinClusEne = g.fMinClusEne;
\r
192 fPairDcaCut = g.fPairDcaCut;
\r
193 fDecayLenCut = g.fDecayLenCut;
\r
194 fImpactCut = g.fImpactCut;
\r
195 fAssocPtCut = g.fAssocPtCut;
\r
196 fMassCut = g.fMassCut;
\r
197 fSdcaCut = g.fSdcaCut;
\r
198 fITSCut = g.fITSCut;
\r
199 fNTagTrkCut = g.fNTagTrkCut;
\r
200 fIPSigCut = g.fIPSigCut;
\r
201 fJetEtaCut = g.fJetEtaCut;
\r
202 fJetPhiMin = g.fJetPhiMin;
\r
203 fJetPhiMax = g.fJetPhiMax;
\r
204 fWriteNtuple = g.fWriteNtuple;
\r
206 fhImpactXY = g.fhImpactXY;
\r
207 fhRefMult = g.fhRefMult;
\r
208 fhRefMult2 = g.fhRefMult2;
\r
210 fh3pOverE = g.fh3pOverE;
\r
211 fh3EOverp = g.fh3EOverp;
\r
212 fh3pOverE2 = g.fh3pOverE2;
\r
213 fh3EOverp2 = g.fh3EOverp2;
\r
214 fh3pOverE3 = g.fh3pOverE3;
\r
215 fh3EOverp3 = g.fh3EOverp3;
\r
216 fh2pOverE = g.fh2pOverE;
\r
217 fh2EOverp = g.fh2EOverp;
\r
218 fh2pOverE2 = g.fh2pOverE2;
\r
219 fh2EOverp2 = g.fh2EOverp2;
\r
221 fh2EledEdx = g.fh2EledEdx;
\r
222 fh2MatchdEdx = g.fh2MatchdEdx;
\r
223 fh2dEtadPhi = g.fh2dEtadPhi;
\r
224 fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;
\r
225 fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;
\r
226 fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;
\r
227 fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;
\r
228 fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;
\r
229 fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;
\r
230 //Photonic electron checks
\r
231 fh1OpeningAngle = g.fh1OpeningAngle;
\r
232 fh1MinvPhoton = g.fh1MinvPhoton;
\r
233 //Reconstructed electrons
\r
234 fhPtElectron = g.fhPtElectron;
\r
235 fhPhiElectron = g.fhPhiElectron;
\r
236 fhEtaElectron = g.fhEtaElectron;
\r
237 fhPtNPE = g.fhPtNPE;
\r
238 fhPhiNPE = g.fhPhiNPE;
\r
239 fhEtaNPE = g.fhEtaNPE;
\r
241 fhPhiPE = g.fhPhiPE;
\r
242 fhEtaPE = g.fhEtaPE;
\r
243 //for comparisons with tracking detectors
\r
244 fhPtHadron = g.fhPtHadron; fhPtNPEleTPC = g.fhPtNPEleTPC;
\r
245 fhPtNPEleTPCTRD = g.fhPtNPEleTPCTRD; fhPtNPEleTTE = g.fhPtNPEleTTE;
\r
246 fhPtNPEleEMCAL = g.fhPtNPEleEMCAL;
\r
248 fhDVMBtagCut1 = g.fhDVMBtagCut1;
\r
249 fhDVMBtagCut2 = g.fhDVMBtagCut2;
\r
250 fhDVMBtagCut3 = g.fhDVMBtagCut3;
\r
251 fhDVMBtagQA1 = g.fhDVMBtagQA1;
\r
252 fhDVMBtagQA2 = g.fhDVMBtagQA2;
\r
253 fhDVMBtagQA3 = g.fhDVMBtagQA3;
\r
254 fhDVMBtagQA4 = g.fhDVMBtagQA4;
\r
255 fhDVMBtagQA5 = g.fhDVMBtagQA5;
\r
257 fhIPSigBtagQA1 = g.fhIPSigBtagQA1;
\r
258 fhIPSigBtagQA2 = g.fhIPSigBtagQA2;
\r
259 fhTagJetPt1x4 = g.fhTagJetPt1x4;
\r
260 fhTagJetPt2x3 = g.fhTagJetPt2x3;
\r
261 fhTagJetPt3x2 = g.fhTagJetPt3x2;
\r
262 fhePlusTagJetPt1x4 = g.fhePlusTagJetPt1x4;
\r
263 fhePlusTagJetPt2x3 = g.fhePlusTagJetPt2x3;
\r
264 fhePlusTagJetPt3x2 = g.fhePlusTagJetPt3x2;
\r
266 fhJetType = g.fhJetType;
\r
267 fhLeadJetType = g.fhLeadJetType;
\r
268 fhBJetXsiFF = g.fhBJetXsiFF;
\r
269 fhBJetPtFF = g.fhBJetPtFF;
\r
270 fhBJetEtaPhi = g.fhBJetEtaPhi;
\r
271 fhNonBJetXsiFF = g.fhNonBJetXsiFF;
\r
272 fhNonBJetPtFF = g.fhNonBJetPtFF;
\r
273 fhNonBJetEtaPhi = g.fhNonBJetEtaPhi;
\r
274 /////////////////////////////////////////////////////////////
\r
275 //Histograms that rely on MC info (not filled for real data)
\r
276 fEleNtuple = g.fEleNtuple;
\r
277 //reco electrons from various sources
\r
278 fhPhiConversion = g.fhPhiConversion;
\r
279 fhEtaConversion = g.fhEtaConversion;
\r
280 //for comparisons with tracking detectors
\r
281 fhPtTrack = g.fhPtTrack;
\r
282 fhPtNPEBHadron = g.fhPtNPEBHadron;
\r
283 //for computing efficiency of B-jet tags
\r
284 fhBJetPt1x4 = g.fhBJetPt1x4; fhBJetPt2x3 = g.fhBJetPt2x3;
\r
285 fhBJetPt3x2 = g.fhBJetPt3x2;
\r
286 fhFakeJetPt1x4 = g.fhFakeJetPt1x4; fhFakeJetPt2x3 = g.fhFakeJetPt2x3;
\r
287 fhFakeJetPt3x2 = g.fhFakeJetPt3x2; fhDVMJet = g.fhDVMJet;
\r
288 //MC rate histograms/ntuple
\r
289 fMCEleNtuple = g.fMCEleNtuple; fhMCBJetElePt = g.fhMCBJetElePt;
\r
290 fhMCBHadronElePt = g.fhMCBHadronElePt;
\r
291 fhPtMCHadron = g.fhPtMCHadron; fhPtMCElectron = g.fhPtMCElectron;
\r
292 fhMCXYConversion = g.fhMCXYConversion;
\r
293 fhMCRadPtConversion = g.fhMCRadPtConversion;
\r
300 //____________________________________________________________________________
\r
301 AliAnaElectron::~AliAnaElectron()
\r
307 //________________________________________________________________________
\r
308 TObjString * AliAnaElectron::GetAnalysisCuts()
\r
310 //Save parameters used for analysis
\r
311 TString parList ; //this will be list of parameters used for this analysis.
\r
314 sprintf(onePar,"--- AliAnaElectron ---\n") ;
\r
316 sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ;
\r
318 sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ;
\r
320 sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ;
\r
322 sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ;
\r
324 sprintf(onePar,"fMinClusEne: %f\n",fMinClusEne) ;
\r
326 sprintf(onePar,"---DVM Btagging\n");
\r
328 sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut);
\r
330 sprintf(onePar,"min ITS-hits: %d\n",fITSCut);
\r
332 sprintf(onePar,"max dR (e,h): %f\n",fDrCut);
\r
334 sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut);
\r
336 sprintf(onePar,"max decaylength: %f\n",fDecayLenCut);
\r
338 sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut);
\r
340 sprintf(onePar,"---IPSig Btagging\n");
\r
342 sprintf(onePar,"min tag track: %d\n",fNTagTrkCut);
\r
344 sprintf(onePar,"min IP significance: %f\n",fIPSigCut);
\r
347 //Get parameters set in base class.
\r
348 parList += GetBaseParametersList() ;
\r
350 //Get parameters set in FiducialCut class (not available yet)
\r
351 //parlist += GetFidCut()->GetFidCutParametersList()
\r
353 return new TObjString(parList) ;
\r
357 //________________________________________________________________________
\r
358 TList * AliAnaElectron::GetCreateOutputObjects()
\r
360 // Create histograms to be saved in output file and
\r
361 // store them in outputContainer
\r
362 TList * outputContainer = new TList() ;
\r
363 outputContainer->SetName("ElectronHistos") ;
\r
365 Int_t nptbins = GetHistoPtBins();
\r
366 Int_t nphibins = GetHistoPhiBins();
\r
367 Int_t netabins = GetHistoEtaBins();
\r
368 Float_t ptmax = GetHistoPtMax();
\r
369 Float_t phimax = GetHistoPhiMax();
\r
370 Float_t etamax = GetHistoEtaMax();
\r
371 Float_t ptmin = GetHistoPtMin();
\r
372 Float_t phimin = GetHistoPhiMin();
\r
373 Float_t etamin = GetHistoEtaMin();
\r
376 fhImpactXY = new TH1F("hImpactXY","Impact parameter for all tracks",200,-10,10.);
\r
377 fhRefMult = new TH1F("hRefMult" ,"refmult QA: " ,5000,0,5000);
\r
378 fhRefMult2 = new TH1F("hRefMult2" ,"refmult2 QA: " ,5000,0,5000);
\r
380 outputContainer->Add(fhImpactXY);
\r
381 outputContainer->Add(fhRefMult);
\r
382 outputContainer->Add(fhRefMult2);
\r
385 fh3pOverE = new TH3F("h3pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);
\r
386 fh3EOverp = new TH3F("h3EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);
\r
387 fh3pOverE2 = new TH3F("h3pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);
\r
388 fh3EOverp2 = new TH3F("h3EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);
\r
389 fh3pOverE3 = new TH3F("h3pOverE_Tpc","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);
\r
390 fh3EOverp3 = new TH3F("h3EOverp_Tpc","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);
\r
391 fh2pOverE = new TH2F("h2pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);
\r
392 fh2EOverp = new TH2F("h2EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );
\r
393 fh2pOverE2 = new TH2F("h2pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);
\r
394 fh2EOverp2 = new TH2F("h2EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );
\r
396 fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi());
\r
397 fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.);
\r
398 fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.);
\r
399 fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
\r
400 fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
\r
401 fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
\r
403 fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
\r
404 fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
\r
405 fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax);
\r
406 fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax);
\r
408 outputContainer->Add(fh3pOverE) ;
\r
409 outputContainer->Add(fh3EOverp) ;
\r
410 outputContainer->Add(fh3pOverE2) ;
\r
411 outputContainer->Add(fh3EOverp2) ;
\r
412 outputContainer->Add(fh3pOverE3) ;
\r
413 outputContainer->Add(fh3EOverp3) ;
\r
414 outputContainer->Add(fh2pOverE) ;
\r
415 outputContainer->Add(fh2EOverp) ;
\r
416 outputContainer->Add(fh2pOverE2) ;
\r
417 outputContainer->Add(fh2EOverp2) ;
\r
418 outputContainer->Add(fh1dR) ;
\r
419 outputContainer->Add(fh2EledEdx) ;
\r
420 outputContainer->Add(fh2MatchdEdx) ;
\r
421 outputContainer->Add(fh2dEtadPhi) ;
\r
422 outputContainer->Add(fh2dEtadPhiMatched) ;
\r
423 outputContainer->Add(fh2dEtadPhiUnmatched) ;
\r
424 outputContainer->Add(fh2TrackPVsClusterE) ;
\r
425 outputContainer->Add(fh2TrackPtVsClusterE) ;
\r
426 outputContainer->Add(fh2TrackPhiVsClusterPhi) ;
\r
427 outputContainer->Add(fh2TrackEtaVsClusterEta) ;
\r
429 //photonic electron checks
\r
430 fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between e+e- pairs",100,0.,TMath::Pi());
\r
431 fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of e+e- pairs",200,0.,2.);
\r
433 outputContainer->Add(fh1OpeningAngle);
\r
434 outputContainer->Add(fh1MinvPhoton);
\r
436 //Reconstructed electrons
\r
437 fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);
\r
438 fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
\r
439 fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
\r
440 fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);
\r
441 fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
\r
442 fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
\r
443 fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);
\r
444 fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
\r
445 fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
\r
447 outputContainer->Add(fhPtElectron) ;
\r
448 outputContainer->Add(fhPhiElectron) ;
\r
449 outputContainer->Add(fhEtaElectron) ;
\r
450 outputContainer->Add(fhPtNPE) ;
\r
451 outputContainer->Add(fhPhiNPE) ;
\r
452 outputContainer->Add(fhEtaNPE) ;
\r
453 outputContainer->Add(fhPtPE) ;
\r
454 outputContainer->Add(fhPhiPE) ;
\r
455 outputContainer->Add(fhEtaPE) ;
\r
457 //These histograms are mixed REAL/MC:
\r
458 //Bins along y-axis are:
\r
459 //0 - unfiltered (filled for both real and MC data)
\r
460 //1 - bottom, 2 - charm, 3 - charm from bottom (MC only)
\r
461 //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown (MC only)
\r
462 //8 - misidentified (MC only)
\r
464 //histograms for comparison to tracking detectors
\r
465 fhPtHadron = new TH2F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
466 fhPtNPEleTPC = new TH2F("hPtNPEleTPC","Non-phot. Electrons identified by TPC w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
467 fhPtNPEleTPCTRD = new TH2F("hPtNPEleTPCTRD","Non-phot. Electrons identified by TPC+TRD w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
468 fhPtNPEleTTE = new TH2F("hPtNPEleTTE","Non-phot. Electrons identified by TPC+TRD+EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
469 fhPtNPEleEMCAL = new TH2F("hPtNPEleEMCAL","Non-phot. Electrons identified by EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
471 outputContainer->Add(fhPtHadron);
\r
472 outputContainer->Add(fhPtNPEleTPC);
\r
473 outputContainer->Add(fhPtNPEleTPCTRD);
\r
474 outputContainer->Add(fhPtNPEleTTE);
\r
475 outputContainer->Add(fhPtNPEleEMCAL);
\r
478 fhDVMBtagCut1 = new TH2F("hdvmbtag_cut1","DVM B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);
\r
479 fhDVMBtagCut2 = new TH2F("hdvmbtag_cut2","DVM B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);
\r
480 fhDVMBtagCut3 = new TH2F("hdvmbtag_cut3","DVM B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);
\r
481 fhDVMBtagQA1 = new TH2F("hdvmbtag_qa1" ,"DVM B-tag QA: pairDCA vs length", 100,0,0.2 ,100,0,1.0);
\r
482 fhDVMBtagQA2 = new TH2F("hdvmbtag_qa2" ,"DVM B-tag QA: signDCA vs mass" , 200,-0.5,0.5 ,100,0,10);
\r
483 fhDVMBtagQA3 = new TH1F("hdvmbtag_qa3" ,"DVM B-tag QA: ITS-Hits electron" ,7,0,7);
\r
484 fhDVMBtagQA4 = new TH1F("hdvmbtag_qa4" ,"DVM B-tag QA: IP d electron" ,200,-3,3);
\r
485 fhDVMBtagQA5 = new TH1F("hdvmbtag_qa5" ,"DVM B-tag QA: IP z electron" ,200,-3,3);
\r
487 outputContainer->Add(fhDVMBtagCut1) ;
\r
488 outputContainer->Add(fhDVMBtagCut2) ;
\r
489 outputContainer->Add(fhDVMBtagCut3) ;
\r
490 outputContainer->Add(fhDVMBtagQA1) ;
\r
491 outputContainer->Add(fhDVMBtagQA2) ;
\r
492 outputContainer->Add(fhDVMBtagQA3) ;
\r
493 outputContainer->Add(fhDVMBtagQA4) ;
\r
494 outputContainer->Add(fhDVMBtagQA5) ;
\r
497 fhIPSigBtagQA1 = new TH1F("hipsigbtag_qa1" ,"IPSig B-tag QA: # tag tracks", 20,0,20);
\r
498 fhIPSigBtagQA2 = new TH1F("hipsigbtag_qa2" ,"IPSig B-tag QA: IP significance", 200,-10.,10.);
\r
499 fhTagJetPt1x4 = new TH1F("hTagJetPt1x4","tagged jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);
\r
500 fhTagJetPt2x3 = new TH1F("hTagJetPt2x3","tagged jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);
\r
501 fhTagJetPt3x2 = new TH1F("hTagJetPt3x2","tagged jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);
\r
502 fhePlusTagJetPt1x4 = new TH1F("hePlusTagJetPt1x4","tagged eJet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);
\r
503 fhePlusTagJetPt2x3 = new TH1F("hePlusTagJetPt2x3","tagged eJet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);
\r
504 fhePlusTagJetPt3x2 = new TH1F("hePlusTagJetPt3x2","tagged eJet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);
\r
506 outputContainer->Add(fhIPSigBtagQA1) ;
\r
507 outputContainer->Add(fhIPSigBtagQA2) ;
\r
508 outputContainer->Add(fhTagJetPt1x4);
\r
509 outputContainer->Add(fhTagJetPt2x3);
\r
510 outputContainer->Add(fhTagJetPt3x2);
\r
511 outputContainer->Add(fhePlusTagJetPt1x4);
\r
512 outputContainer->Add(fhePlusTagJetPt2x3);
\r
513 outputContainer->Add(fhePlusTagJetPt3x2);
\r
516 fhJetType = new TH2F("hJetType","# jets passing each tag method vs jet pt",15,0,15,300,0.,300.);
\r
517 fhLeadJetType = new TH2F("hLeadJetType","# leading jets passing each tag method vs jet pt",15,0,15,300,0.,300.);
\r
518 fhBJetXsiFF = new TH2F("hBJetXsiFF","B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);
\r
519 fhBJetPtFF = new TH2F("hBJetPtFF","B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);
\r
520 fhBJetEtaPhi = new TH2F("hBJetEtaPhi","B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);
\r
521 fhNonBJetXsiFF = new TH2F("hNonBJetXsiFF","Non B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);
\r
522 fhNonBJetPtFF = new TH2F("hNonBJetPtFF","Non B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);
\r
523 fhNonBJetEtaPhi = new TH2F("hNonBJetEtaPhi","Non B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);
\r
525 outputContainer->Add(fhJetType);
\r
526 outputContainer->Add(fhLeadJetType);
\r
527 outputContainer->Add(fhBJetXsiFF);
\r
528 outputContainer->Add(fhBJetPtFF);
\r
529 outputContainer->Add(fhBJetEtaPhi);
\r
530 outputContainer->Add(fhNonBJetXsiFF);
\r
531 outputContainer->Add(fhNonBJetPtFF);
\r
532 outputContainer->Add(fhNonBJetEtaPhi);
\r
534 //Histograms that use MC information
\r
537 //electron ntuple for further analysis
\r
539 fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");
\r
540 outputContainer->Add(fEleNtuple) ;
\r
543 //electrons from various MC sources
\r
544 fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
\r
545 fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. pT",nptbins,ptmin,ptmax,netabins,etamin,etamax);
\r
547 outputContainer->Add(fhPhiConversion);
\r
548 outputContainer->Add(fhEtaConversion);
\r
550 //Bins along y-axis are: 0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,
\r
551 //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown, 8 - misidentified
\r
553 //histograms for comparison to tracking detectors
\r
554 fhPtTrack = new TH2F("hPtTrack","Track w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
555 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
557 outputContainer->Add(fhPtTrack);
\r
558 outputContainer->Add(fhPtNPEBHadron);
\r
560 //for computing efficiency of IPSig tag
\r
561 fhBJetPt1x4 = new TH1F("hBJetPt1x4","tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);
\r
562 fhBJetPt2x3 = new TH1F("hBJetPt2x3","tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);
\r
563 fhBJetPt3x2 = new TH1F("hBJetPt3x2","tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);
\r
564 fhFakeJetPt1x4 = new TH1F("hFakeJetPt1x4","fake tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);
\r
565 fhFakeJetPt2x3 = new TH1F("hFakeJetPt2x3","fake tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);
\r
566 fhFakeJetPt3x2 = new TH1F("hFakeJetPt3x2","fake tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);
\r
567 fhDVMJet = new TH2F("hDVM_algo","# DVM jets passing vs Mc-Bjet",10,0,10,300,0.,300.);
\r
569 outputContainer->Add(fhBJetPt1x4);
\r
570 outputContainer->Add(fhBJetPt2x3);
\r
571 outputContainer->Add(fhBJetPt3x2);
\r
572 outputContainer->Add(fhFakeJetPt1x4);
\r
573 outputContainer->Add(fhFakeJetPt2x3);
\r
574 outputContainer->Add(fhFakeJetPt3x2);
\r
575 outputContainer->Add(fhDVMJet);
\r
577 //MC Only histograms
\r
579 //MC ele ntuple for further analysis
\r
581 fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");
\r
582 outputContainer->Add(fMCEleNtuple) ;
\r
585 fhMCBJetElePt = new TH2F("hMCBJetElePt","MC B-jet pT vs. electron pT",300,0.,300.,300,0.,300.);
\r
586 fhMCBHadronElePt = new TH2F("hMCBHadronElePt","MC B-hadron pT vs. electron pT",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
\r
587 fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);
\r
589 //Bins along y-axis are: 0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,
\r
590 //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown
\r
591 fhPtMCElectron = new TH2F("hPtMCElectron","MC electrons from various sources w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);
\r
593 fhMCXYConversion = new TH2F("hMCXYConversion","XvsY of conversion electrons",400,-400.,400.,400,-400.,400.);
\r
594 fhMCRadPtConversion = new TH2F("hMCRadPtConversion","Radius vs pT of conversion electrons",200,0.,400.,nptbins,ptmin,ptmax);
\r
596 outputContainer->Add(fhMCBJetElePt);
\r
597 outputContainer->Add(fhMCBHadronElePt);
\r
598 outputContainer->Add(fhPtMCHadron);
\r
599 outputContainer->Add(fhPtMCElectron);
\r
600 outputContainer->Add(fhMCXYConversion);
\r
601 outputContainer->Add(fhMCRadPtConversion);
\r
606 return outputContainer ;
\r
610 //____________________________________________________________________________
\r
611 void AliAnaElectron::Init()
\r
614 //do some initialization
\r
615 if(fCalorimeter == "PHOS") {
\r
616 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");
\r
617 fCalorimeter = "EMCAL";
\r
619 if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
\r
620 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");
\r
627 //____________________________________________________________________________
\r
628 void AliAnaElectron::InitParameters()
\r
631 //Initialize the parameters of the analysis.
\r
632 SetOutputAODClassName("AliAODPWG4Particle");
\r
633 SetOutputAODName("PWG4Particle");
\r
635 AddToHistogramsName("AnaElectron_");
\r
637 fCalorimeter = "EMCAL" ;
\r
640 fResidualCut = 0.02;
\r
644 fPairDcaCut = 0.02;
\r
645 fDecayLenCut = 1.0;
\r
654 //Jet fiducial cuts
\r
660 //__________________________________________________________________
\r
661 void AliAnaElectron::MakeAnalysisFillAOD()
\r
664 // Do analysis and fill aods with electron candidates
\r
665 // These AODs will be used to do subsequent histogram filling
\r
667 // Also fill some QA histograms
\r
670 Double_t bfield = 0.;
\r
671 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
\r
673 //Select the calorimeter of the electron
\r
674 if(fCalorimeter != "EMCAL") {
\r
675 printf("This class not yet implemented for PHOS\n");
\r
679 TObjArray *cl = GetAODEMCAL();
\r
681 ////////////////////////////////////////////////
\r
682 //Start from tracks and get associated clusters
\r
683 ////////////////////////////////////////////////
\r
684 if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
\r
685 Int_t ntracks = GetAODCTS()->GetEntriesFast();
\r
686 Int_t refmult = 0; Int_t refmult2 = 0;
\r
688 printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
\r
690 //Unfortunately, AliAODTracks don't have associated EMCAL clusters.
\r
691 //we have to redo track-matching, I guess
\r
692 Int_t iCluster = -999;
\r
693 Int_t bt = 0; //counter for event b-tags
\r
695 for (Int_t itrk = 0; itrk < ntracks; itrk++) {////////////// track loop
\r
696 iCluster = -999; //start with no match
\r
697 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;
\r
698 if (TMath::Abs(track->Eta())< 0.5) refmult++;
\r
699 Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};
\r
700 Bool_t dcaOkay = GetDCA(track,imp,cov); //homegrown dca calculation until AOD is fixed
\r
701 if(!dcaOkay) printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d. Skipping it...\n",itrk);
\r
702 if(TMath::Abs(track->Eta())< 0.5 && TMath::Abs(imp[0])<1.0 && TMath::Abs(imp[1])<1.0) refmult2++;
\r
703 fhImpactXY->Fill(imp[0]);
\r
706 //AliESDtrack esdTrack(track);
\r
707 //Double_t tpcpid[AliPID::kSPECIES];
\r
708 //esdTrack.GetTPCpid(tpcpid);
\r
709 //Double_t eProb = tpcpid[AliPID::kElectron];
\r
710 //if(eProb > 0) printf("<%d> ESD eProb = %2.2f\n",itrk,eProb);
\r
712 AliAODPid* pid = (AliAODPid*) track->GetDetPid();
\r
714 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);
\r
717 Double_t emcpos[3];
\r
718 pid->GetEMCALPosition(emcpos);
\r
719 Double_t emcmom[3];
\r
720 pid->GetEMCALMomentum(emcmom);
\r
722 TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);
\r
723 TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);
\r
724 Double_t tphi = pos.Phi();
\r
725 Double_t teta = pos.Eta();
\r
726 Double_t tmom = mom.Mag();
\r
728 //TLorentzVector mom2(mom,0.);
\r
729 Bool_t in = kFALSE;
\r
730 if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&
\r
731 mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE;
\r
732 //Also check the track
\r
733 if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. &&
\r
734 track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE;
\r
735 ////////////////////////////
\r
736 //THIS HAS A MEM LEAK JLK 24-Oct-09
\r
737 //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
\r
738 ///////////////////////////
\r
739 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Track(Extrap) pt %2.2f(%2.2f), phi %2.2f(%2.2f), eta %2.2f(%2.2f) in fiducial cut %d\n",track->Pt(), mom.Pt(), track->Phi(), mom.Phi(), track->Eta(),mom.Eta(), in);
\r
741 if(mom.Pt() > GetMinPt() && in) {
\r
743 Double_t dEdx = pid->GetTPCsignal();
\r
745 Int_t ntot = cl->GetEntriesFast();
\r
746 Double_t res = 999.;
\r
747 Double_t pOverE = -999.;
\r
749 Int_t pidProb = track->GetMostProbablePID();
\r
750 Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;
\r
751 Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;
\r
752 Bool_t trkChgHad = kFALSE; if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) trkChgHad = kTRUE;
\r
756 //Check against V0 for conversion, only if it is flagged as electron
\r
757 Bool_t photonic = kFALSE;
\r
758 if(tpcEle || trkEle) photonic = PhotonicV0(itrk);
\r
759 if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),0); //0 = no MC info
\r
760 if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),0); //0 = no MC info
\r
762 if(trkChgHad) fhPtHadron->Fill(track->Pt(),0); //0 = no MC info
\r
764 //Input from second AOD?
\r
766 //if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;
\r
767 tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);
\r
769 if(trkChgHad) fhPtHadron->Fill(track->Pt(),GetMCSource(tmctag));
\r
770 if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),GetMCSource(tmctag));
\r
771 if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),GetMCSource(tmctag));
\r
772 fhPtTrack->Fill(track->Pt(),GetMCSource(tmctag));
\r
775 Bool_t emcEle = kFALSE;
\r
776 //For tracks in EMCAL acceptance, pair them with all clusters
\r
777 //and fill the dEta vs dPhi for these pairs:
\r
779 Double_t minR = 99;
\r
780 Double_t minPe =-1;
\r
781 Double_t minEp =-1;
\r
782 Double_t minMult = -1;
\r
783 Double_t minPt = -1;
\r
785 for(Int_t iclus = 0; iclus < ntot; iclus++) {
\r
786 AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));
\r
787 if(!clus) continue;
\r
789 //As of 11-Oct-2009
\r
790 //only select "good" clusters
\r
791 if (clus->GetNCells() < 2 ) continue;
\r
792 if (clus->GetNCells() > 30 ) continue;
\r
793 if (clus->E() < fMinClusEne ) continue;
\r
794 if (clus->GetDispersion() > 1 ) continue;
\r
795 if (clus->GetM20() > 0.4 ) continue;
\r
796 if (clus->GetM02() > 0.4 ) continue;
\r
797 if (clus->GetM20() < 0.03 ) continue;
\r
798 if (clus->GetM02() < 0.03 ) continue;
\r
801 clus->GetPosition(x);
\r
802 TVector3 cluspos(x[0],x[1],x[2]);
\r
803 Double_t deta = teta - cluspos.Eta();
\r
804 Double_t dphi = tphi - cluspos.Phi();
\r
805 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
\r
806 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
\r
807 fh2dEtadPhi->Fill(deta,dphi);
\r
808 fh2TrackPVsClusterE->Fill(clus->E(),track->P());
\r
809 fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());
\r
810 fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());
\r
811 fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());
\r
813 res = sqrt(dphi*dphi + deta*deta);
\r
816 /////////////////////////////////
\r
817 //Perform electron cut analysis//
\r
818 /////////////////////////////////
\r
820 if(res < fResidualCut) {
\r
821 fh2dEtadPhiMatched->Fill(deta,dphi);
\r
822 fh2MatchdEdx->Fill(track->P(),dEdx);
\r
825 Double_t energy = clus->E();
\r
826 if(energy > 0) pOverE = tmom/energy;
\r
831 minEp = energy/tmom;
\r
832 minMult = clus->GetNCells() ;
\r
833 minPt = track->Pt();
\r
836 Int_t cmctag = -1;
\r
838 //Do you want the cluster or the track label?
\r
840 //if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;
\r
841 cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(),GetReader(),input);
\r
845 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
849 fh2dEtadPhiUnmatched->Fill(deta,dphi);
\r
851 }//calo cluster loop
\r
853 fh3pOverE->Fill(minPt,minPe ,minMult);
\r
854 fh3EOverp->Fill(minPt,minEp ,minMult);
\r
856 fh3pOverE2->Fill(minPt,minPe ,minMult);
\r
857 fh3EOverp2->Fill(minPt,minEp ,minMult);
\r
860 fh3pOverE3->Fill(minPt,minPe ,minMult);
\r
861 fh3EOverp3->Fill(minPt,minEp ,minMult);
\r
864 if (tmctag>-1 && GetMCSource(tmctag)<8 ) {
\r
865 fh2pOverE->Fill(minPt,minPe );
\r
866 fh2EOverp->Fill(minPt,minEp );
\r
868 fh2pOverE2->Fill(minPt,minPe );
\r
869 fh2EOverp2->Fill(minPt,minEp );
\r
873 //////////////////////////////
\r
874 //Electron cuts happen here!//
\r
875 //////////////////////////////
\r
876 if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;
\r
878 ///////////////////////////
\r
879 //Fill AOD with electrons//
\r
880 ///////////////////////////
\r
881 //Take all emcal electrons, but the others only if pT < 10 GeV
\r
882 if(emcEle || ( (tpcEle || trkEle) && track->Pt() < 10.) ) {
\r
885 if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");
\r
886 Int_t dvmbtag = GetDVMBtag(track); bt += dvmbtag;
\r
888 fh2EledEdx->Fill(track->P(),dEdx);
\r
890 Double_t eMass = 0.511/1000; //mass in GeV
\r
891 Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);
\r
892 AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);
\r
893 tr.SetLabel(track->GetLabel());
\r
894 tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters
\r
895 tr.SetTrackLabel(itrk,-1); //sets the indices of the original tracks
\r
896 if(emcEle) {//PID determined by EMCAL
\r
897 tr.SetDetector(fCalorimeter);
\r
899 tr.SetDetector("CTS"); //PID determined by CTS
\r
902 //if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);
\r
903 //Make this preserve sign of particle
\r
904 if(track->Charge() < 0) tr.SetPdg(11); //electron is 11
\r
905 else tr.SetPdg(-11); //positron is -11
\r
907 if(dvmbtag > 0) tr.SetBTagBit(btag,tr.kDVMTag0);
\r
908 if(dvmbtag > 1) tr.SetBTagBit(btag,tr.kDVMTag1);
\r
909 if(dvmbtag > 2) tr.SetBTagBit(btag,tr.kDVMTag2);
\r
912 //Play with the MC stack if available
\r
913 //Check origin of the candidates
\r
916 //FIXME: Need to re-think this for track-oriented analysis
\r
917 //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?
\r
918 tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));
\r
920 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());
\r
921 }//Work with stack also
\r
923 AddAODParticle(tr);
\r
925 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());
\r
927 }//pt, fiducial selection
\r
931 fhRefMult->Fill(refmult);
\r
932 fhRefMult2->Fill(refmult2);
\r
934 if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");
\r
935 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs \n");
\r
939 //__________________________________________________________________
\r
940 void AliAnaElectron::MakeAnalysisFillHistograms()
\r
942 //Do analysis and fill histograms
\r
944 AliStack * stack = 0x0;
\r
945 TParticle * primary = 0x0;
\r
946 AliAODMCParticle * aodprimary = 0x0;
\r
948 Int_t ph1 = 0; //photonic 1 count
\r
949 Int_t ph2 = 0; //photonic 2 count
\r
950 Int_t phB = 0; //both count
\r
953 if(GetReader()->ReadStack()){
\r
954 stack = GetMCStack() ;
\r
956 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");
\r
961 ////////////////////////////////////
\r
962 //Loop over jets and check for b-tag
\r
963 ////////////////////////////////////
\r
964 Int_t njets = (GetReader()->GetOutputEvent())->GetNJets();
\r
966 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets. Performing b-jet tag analysis\n",njets);
\r
968 for(Int_t ijet = 0; ijet < njets ; ijet++) {
\r
969 AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;
\r
970 //Only consider jets with pt > 10 GeV (the rest have to be junk)
\r
971 //printf("AODJet<%d> pt = %2.2f\n",ijet,jet->Pt());
\r
972 if(jet->Pt() < 10.) continue;
\r
974 if(GetDebug() > 3) {
\r
975 printf("AliAODJet ijet = %d\n",ijet);
\r
978 //Skip jets not inside a smaller fiducial volume to ensure that
\r
979 //they are completely contained in the EMCAL
\r
980 if(TMath::Abs(jet->Eta()) > fJetEtaCut) continue;
\r
981 if(jet->Phi() < fJetPhiMin || jet->Phi() > fJetPhiMax) continue;
\r
983 //To "tag" the jet, we will look for it to pass our various criteria
\r
984 //For e jet tag, we just look to see which ones have NPEs
\r
985 //For DVM jet tag, we will look for DVM electrons
\r
986 //For IPSig, we compute the IPSig for all tracks and if the
\r
987 //number passing is above the cut, it passes
\r
988 Bool_t leadJet = kFALSE;
\r
989 if (ijet==0) leadJet= kTRUE;
\r
990 Bool_t eJet = kFALSE;
\r
991 Bool_t eJet2 = kFALSE; //electron triggered
\r
992 Bool_t hadJet = kFALSE; //hadron triggered
\r
993 Bool_t dvmJet = kFALSE;
\r
994 Bool_t ipsigJet = kFALSE;
\r
995 TRefArray* rt = jet->GetRefTracks();
\r
996 Int_t ntrk = rt->GetEntries();
\r
997 Int_t trackCounter[4] = {0,0,0,0}; //for ipsig
\r
998 for(Int_t itrk = 0; itrk < ntrk; itrk++) {
\r
999 AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);
\r
1000 if( GetIPSignificance(jetTrack, jet->Phi()) > fIPSigCut) trackCounter[0]++;
\r
1001 if( GetIPSignificance(jetTrack, jet->Phi()) > 4.) trackCounter[1]++;
\r
1002 if( GetIPSignificance(jetTrack, jet->Phi()) > 3.) trackCounter[2]++;
\r
1003 if( GetIPSignificance(jetTrack, jet->Phi()) > 2.) trackCounter[3]++;
\r
1004 Bool_t isNPE = CheckTrack(jetTrack,"NPE");
\r
1005 if(isNPE) eJet = kTRUE;
\r
1006 if ( isNPE && jetTrack->Pt()>10.0 ) eJet2 =kTRUE;
\r
1007 if (!isNPE && jetTrack->Pt()>10.0) hadJet =kTRUE;
\r
1008 Bool_t isDVM = CheckTrack(jetTrack,"DVM");
\r
1009 if(isDVM) dvmJet = kTRUE;
\r
1011 fhIPSigBtagQA1->Fill(trackCounter[0]);
\r
1012 if(trackCounter[1]>0) fhTagJetPt1x4->Fill(jet->Pt());
\r
1013 if(trackCounter[2]>1) fhTagJetPt2x3->Fill(jet->Pt());
\r
1014 if(trackCounter[3]>2) fhTagJetPt3x2->Fill(jet->Pt());
\r
1016 if(trackCounter[1]>0 && eJet) fhePlusTagJetPt1x4->Fill(jet->Pt());
\r
1017 if(trackCounter[2]>1 && eJet) fhePlusTagJetPt2x3->Fill(jet->Pt());
\r
1018 if(trackCounter[3]>2 && eJet) fhePlusTagJetPt3x2->Fill(jet->Pt());
\r
1020 if(trackCounter[0] > fNTagTrkCut) ipsigJet = kTRUE;
\r
1023 //determine tagging efficiency & mis-tagging rate
\r
1024 //using b-quarks from stack
\r
1025 Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());
\r
1026 Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());
\r
1027 if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");
\r
1028 if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");
\r
1029 if (dvmJet && GetDebug() > 0) printf("== found DVM jet==\n");
\r
1031 if(isTrueBjet && dvmJet) fhDVMJet->Fill(0.,jet->Pt()); // good tagged
\r
1032 if(isTrueBjet && !dvmJet) fhDVMJet->Fill(1.,jet->Pt()); // missed tagged
\r
1033 if(!isTrueBjet && dvmJet) fhDVMJet->Fill(2.,jet->Pt()); // fake tagged
\r
1034 if(!isTrueBjet && !dvmJet) fhDVMJet->Fill(3.,jet->Pt()); // others
\r
1035 if(isTrueDjet && !isTrueBjet && dvmJet) fhDVMJet->Fill(4.,jet->Pt()); // charm-tagged
\r
1036 if(isTrueDjet && !isTrueBjet && !dvmJet) fhDVMJet->Fill(5.,jet->Pt()); // charm -not tagged
\r
1037 if(!(isTrueDjet||isTrueBjet ) && dvmJet) fhDVMJet->Fill(6.,jet->Pt()); // light flavor -tagged
\r
1038 if(!(isTrueDjet||isTrueBjet ) && !dvmJet) fhDVMJet->Fill(7.,jet->Pt()); // light flavor -not tagged
\r
1039 if(isTrueBjet && eJet && dvmJet) fhDVMJet->Fill(8.,jet->Pt()); // bjet with electron
\r
1040 if(isTrueBjet && !eJet && dvmJet) fhDVMJet->Fill(9.,jet->Pt()); // needs more thought
\r
1043 if(trackCounter[1]>0) fhBJetPt1x4->Fill(jet->Pt());
\r
1044 if(trackCounter[2]>1) fhBJetPt2x3->Fill(jet->Pt());
\r
1045 if(trackCounter[3]>2) fhBJetPt3x2->Fill(jet->Pt());
\r
1047 if(trackCounter[1]>0) fhFakeJetPt1x4->Fill(jet->Pt());
\r
1048 if(trackCounter[2]>1) fhFakeJetPt2x3->Fill(jet->Pt());
\r
1049 if(trackCounter[3]>2) fhFakeJetPt3x2->Fill(jet->Pt());
\r
1053 //Fill bjet histograms here
\r
1054 if(!(eJet || ipsigJet || dvmJet)) fhJetType->Fill(0.,jet->Pt()); //none
\r
1055 if(eJet && !(ipsigJet || dvmJet)) fhJetType->Fill(1.,jet->Pt()); //only ejet
\r
1056 if(dvmJet && !(eJet || ipsigJet)) fhJetType->Fill(2.,jet->Pt()); //only dvm
\r
1057 if(ipsigJet && !(eJet || dvmJet)) fhJetType->Fill(3.,jet->Pt()); //only ipsig
\r
1058 if(eJet && dvmJet && !ipsigJet) fhJetType->Fill(4.,jet->Pt()); //ejet & dvm
\r
1059 if(eJet && ipsigJet && !dvmJet) fhJetType->Fill(5.,jet->Pt()); //ejet & ipsig
\r
1060 if(dvmJet && ipsigJet && !eJet) fhJetType->Fill(6.,jet->Pt()); //dvm & ipsig
\r
1061 if(dvmJet && ipsigJet && eJet) fhJetType->Fill(7.,jet->Pt()); //all
\r
1062 if(dvmJet || ipsigJet || eJet) fhJetType->Fill(8.,jet->Pt()); //any of them
\r
1063 if(eJet ) fhJetType->Fill(9.,jet->Pt()); //any of them
\r
1064 if(dvmJet) fhJetType->Fill(10.,jet->Pt()); //any of them
\r
1065 if(eJet2 ) fhJetType->Fill(11.,jet->Pt()); //any of them
\r
1066 if(hadJet) fhJetType->Fill(12.,jet->Pt()); //any of them
\r
1068 if(eJet || ipsigJet || dvmJet) fhBJetEtaPhi->Fill(jet->Eta(),jet->Phi());
\r
1069 else fhNonBJetEtaPhi->Fill(jet->Eta(),jet->Phi());
\r
1073 fhLeadJetType->Fill(0.,jet->Pt()); //all
\r
1074 if(eJet ) fhLeadJetType->Fill(1.,jet->Pt());
\r
1075 if(eJet2 ) fhLeadJetType->Fill(2.,jet->Pt());
\r
1076 if(hadJet ) fhLeadJetType->Fill(3.,jet->Pt());
\r
1077 if(eJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(4.,jet->Pt());
\r
1078 if(eJet2 && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(5.,jet->Pt());
\r
1079 if(hadJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(6.,jet->Pt());
\r
1082 for(Int_t itrk = 0; itrk < ntrk; itrk++) {
\r
1083 AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);
\r
1084 Double_t xsi = TMath::Log(jet->Pt()/jetTrack->Pt());
\r
1085 if(eJet || ipsigJet || dvmJet) {
\r
1086 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms - We have a bjet!\n");
\r
1087 fhBJetXsiFF->Fill(xsi,jet->Pt());
\r
1088 fhBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());
\r
1090 //Fill non-bjet histograms here
\r
1091 fhNonBJetXsiFF->Fill(xsi,jet->Pt());
\r
1092 fhNonBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());
\r
1099 //////////////////////////////
\r
1100 //Loop on stored AOD electrons
\r
1101 //////////////////////////////
\r
1102 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
\r
1103 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
\r
1105 for(Int_t iaod = 0; iaod < naod ; iaod++){
\r
1106 AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
\r
1107 Int_t pdg = ele->GetPdg();
\r
1109 if(GetDebug() > 3)
\r
1110 printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;
\r
1112 if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue;
\r
1114 if(GetDebug() > 1)
\r
1115 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;
\r
1117 //MC tag of this electron
\r
1118 Int_t mctag = ele->GetTag();
\r
1120 //Filter for photonic electrons based on opening angle and Minv
\r
1121 //cuts, also fill histograms
\r
1122 Bool_t photonic = kFALSE;
\r
1123 Bool_t photonic1 = kFALSE;
\r
1124 photonic1 = PhotonicPrim(ele); //check against primaries
\r
1125 if(photonic1) ph1++;
\r
1126 Bool_t photonic2 = kFALSE;
\r
1127 photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s
\r
1128 if(photonic2) ph2++;
\r
1129 if(photonic1 && photonic2) phB++;
\r
1130 if(photonic1 || photonic2) photonic = kTRUE;
\r
1132 //Fill electron histograms
\r
1133 Float_t ptele = ele->Pt();
\r
1134 Float_t phiele = ele->Phi();
\r
1135 Float_t etaele = ele->Eta();
\r
1137 //"Best reconstructed electron spectrum" = EMCAL or tracking
\r
1138 //detectors say it is an electron and it does not form a V0
\r
1139 //with Minv near a relevant resonance
\r
1141 fhPtNPEleTTE->Fill(ptele,0); //0 = no MC info
\r
1142 if(ele->GetDetector() == fCalorimeter) fhPtNPEleEMCAL->Fill(ptele,0);
\r
1144 fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag));
\r
1145 if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag));
\r
1146 if(GetMCSource(mctag) == 1) { //it's a bottom electron, now
\r
1147 //get the parent's pt
\r
1148 Double_t ptbHadron = GetBParentPt(ele->GetLabel());
\r
1149 fhPtNPEBHadron->Fill(ptele,ptbHadron);
\r
1154 //kept for historical reasons?
\r
1155 fhPtElectron ->Fill(ptele);
\r
1156 fhPhiElectron ->Fill(ptele,phiele);
\r
1157 fhEtaElectron ->Fill(ptele,etaele);
\r
1160 fhPtPE->Fill(ptele);
\r
1161 fhPhiPE->Fill(ptele,phiele);
\r
1162 fhEtaPE->Fill(ptele,etaele);
\r
1164 fhPtNPE->Fill(ptele);
\r
1165 fhPhiNPE->Fill(ptele,phiele);
\r
1166 fhEtaNPE->Fill(ptele,etaele);
\r
1170 if(GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCConversion)){
\r
1171 fhPhiConversion ->Fill(ptele,phiele);
\r
1172 fhEtaConversion ->Fill(ptele,etaele);
\r
1174 }//Histograms with MC
\r
1178 ////////////////////////////////////////////////////////
\r
1179 //Fill histograms of pure MC kinematics from the stack//
\r
1180 ////////////////////////////////////////////////////////
\r
1184 TVector3 bjetVect[4];
\r
1185 Int_t nPythiaGenJets = 0;
\r
1186 AliGenPythiaEventHeader* pythiaGenHeader = (AliGenPythiaEventHeader*)GetReader()->GetGenEventHeader();
\r
1187 if(pythiaGenHeader){
\r
1188 //Get Jets from MC header
\r
1189 nPythiaGenJets = pythiaGenHeader->NTriggerJets();
\r
1191 for(int ip = 0;ip < nPythiaGenJets;++ip){
\r
1192 if (iCount>3) break;
\r
1194 pythiaGenHeader->TriggerJet(ip,p);
\r
1195 TVector3 tempVect(p[0],p[1],p[2]);
\r
1196 if ( TMath::Abs(tempVect.Eta())>fJetEtaCut || tempVect.Phi() < fJetPhiMin || tempVect.Phi() > fJetPhiMax) continue;
\r
1197 //Only store it if it has a b-quark within dR < 0.2 of jet axis ?
\r
1198 if(IsMcBJet(tempVect.Eta(),tempVect.Phi())) {
\r
1199 bjetVect[iCount].SetXYZ(p[0], p[1], p[2]);
\r
1205 Int_t nPart = GetNumAODMCParticles();
\r
1206 if(GetReader()->ReadStack()) nPart = stack->GetNtrack();
\r
1208 for(Int_t ipart = 0; ipart < nPart; ipart++) {
\r
1210 //All the variables we want from MC particles
\r
1211 Double_t px = 0.; Double_t py = 0.; Double_t pz = 0.; Double_t e = 0.;
\r
1212 Double_t vx = -999.; Double_t vy = -999.; Double_t vz = -999.; Double_t vt = -999.;
\r
1213 Int_t pdg = 0; Int_t mpdg = 0; Double_t mpt = 0.;
\r
1215 if(GetReader()->ReadStack()) {
\r
1216 primary = stack->Particle(ipart);
\r
1217 pdg = primary->GetPdgCode();
\r
1218 px = primary->Px(); py = primary->Py(); pz = primary->Pz(); e = primary->Energy();
\r
1219 vx = primary->Vx(); vy = primary->Vy(); vz = primary->Vz(); vt = primary->T();
\r
1220 if(primary->GetMother(0)>=0) {
\r
1221 TParticle *parent = stack->Particle(primary->GetMother(0));
\r
1223 mpdg = parent->GetPdgCode();
\r
1224 mpt = parent->Pt();
\r
1227 } else if(GetReader()->ReadAODMCParticles()) {
\r
1228 aodprimary = (AliAODMCParticle*)GetMCParticle(ipart);
\r
1229 pdg = aodprimary->GetPdgCode();
\r
1230 px = aodprimary->Px(); py = aodprimary->Py(); pz = aodprimary->Pz(); e = aodprimary->E();
\r
1231 vx = aodprimary->Xv(); vy = aodprimary->Yv(); vz = aodprimary->Zv(); vt = aodprimary->T();
\r
1232 Int_t parentId = aodprimary->GetMother();
\r
1234 AliAODMCParticle *parent = (AliAODMCParticle*)GetMCParticle(parentId);
\r
1236 mpdg = parent->GetPdgCode();
\r
1237 mpt = parent->Pt();
\r
1242 TLorentzVector mom(px,py,pz,e);
\r
1243 TLorentzVector pos(vx,vy,vz,vt);
\r
1244 Bool_t in = kFALSE;
\r
1245 if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&
\r
1246 mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE;
\r
1247 /////////////////////////////////
\r
1248 //THIS HAS A MEM LEAK JLK 24-Oct-09
\r
1249 //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter);
\r
1250 ////////////////////////////////
\r
1251 if(mom.Pt() < GetMinPt()) continue;
\r
1254 if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 2212)
\r
1255 fhPtMCHadron->Fill(mom.Pt());
\r
1257 //we only care about electrons
\r
1258 if(TMath::Abs(pdg) != 11) continue;
\r
1259 //we only want TRACKABLE electrons (TPC 85-250cm)
\r
1260 if(pos.Rho() > 200.) continue;
\r
1261 //Ignore low pt electrons
\r
1262 if(mom.Pt() < 0.2) continue;
\r
1264 //find out what the ancestry of this electron is
\r
1267 mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);
\r
1269 if(GetMCSource(mctag)==1) { //bottom electron
\r
1270 //See if it is within dR < 0.4 of a bjet
\r
1271 for(Int_t ij = 0; ij < nPythiaGenJets; ij++) {
\r
1272 Double_t deta = primary->Eta() - bjetVect[ij].Eta();
\r
1273 Double_t dphi = primary->Phi() - bjetVect[ij].Phi();
\r
1274 Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi);
\r
1276 fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt());
\r
1281 if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) ||
\r
1282 (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )
\r
1284 fhMCBHadronElePt->Fill(mom.Pt(), mpt);
\r
1286 //CHECK THAT THIS IS CORRECTLY FILLED - SHOULD WE USE MCSOURCE HERE?
\r
1287 fhPtMCElectron->Fill(mom.Pt(),0); //0 = unfiltered
\r
1288 fhPtMCElectron->Fill(mom.Pt(),GetMCSource(mctag));
\r
1290 if(GetMCSource(mctag) == 4) {//conversion
\r
1291 fhMCXYConversion->Fill(vx,vy);
\r
1292 fhMCRadPtConversion->Fill(TMath::Sqrt(vx*vx+vy*vy),mom.Pt());
\r
1296 if(fWriteNtuple) {
\r
1297 fMCEleNtuple->Fill(mctag,mom.Pt(),mom.Phi(),mom.Eta(),vx,vy,vz);
\r
1302 //if(GetDebug() > 0)
\r
1303 printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB);
\r
1306 //__________________________________________________________________
\r
1307 Int_t AliAnaElectron::GetDVMBtag(AliAODTrack * tr )
\r
1309 //This method uses the Displaced Vertex between electron-hadron
\r
1310 //pairs and the primary vertex to determine whether an electron is
\r
1311 //likely from a B hadron.
\r
1314 for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;
\r
1316 fhDVMBtagQA3->Fill(ncls1);
\r
1317 if (ncls1 < fITSCut) return 0;
\r
1319 Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};
\r
1320 Bool_t dcaOkay = GetDCA(tr,imp,cov); //homegrown dca calculation until AOD is fixed
\r
1322 printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d",tr->GetID());
\r
1326 fhDVMBtagQA4->Fill(imp[0]);
\r
1327 if (TMath::Abs(imp[0]) > fImpactCut ) return 0;
\r
1328 fhDVMBtagQA5->Fill(imp[1]);
\r
1329 if (TMath::Abs(imp[1]) > fImpactCut ) return 0;
\r
1335 for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {
\r
1337 AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);
\r
1338 Int_t id1 = tr->GetID();
\r
1339 Int_t id2 = track2->GetID();
\r
1340 if(id1 == id2) continue;
\r
1343 for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
\r
1344 if (ncls2 < fITSCut) continue;
\r
1346 if(track2->Pt() < fAssocPtCut) continue;
\r
1348 Double_t dphi = tr->Phi() - track2->Phi();
\r
1349 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
\r
1350 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
\r
1351 Double_t deta = tr->Eta() - track2->Eta();
\r
1352 Double_t dr = sqrt(deta*deta + dphi*dphi);
\r
1354 if(dr > fDrCut) continue;
\r
1356 Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);
\r
1357 if (sDca1 > fSdcaCut) nvtx1++;
\r
1358 Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);
\r
1359 if (sDca2 > fSdcaCut) nvtx2++;
\r
1360 Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);
\r
1361 if (sDca3 > fSdcaCut) nvtx3++;
\r
1363 } //loop over hadrons
\r
1365 if(GetDebug() > 0) {
\r
1366 if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);
\r
1367 if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);
\r
1368 if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);
\r
1371 //fill QA histograms
\r
1372 fhDVMBtagCut1->Fill(nvtx1,tr->Pt());
\r
1373 fhDVMBtagCut2->Fill(nvtx2,tr->Pt());
\r
1374 fhDVMBtagCut3->Fill(nvtx3,tr->Pt());
\r
1380 //__________________________________________________________________
\r
1381 Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)
\r
1383 //Compute the signed dca between two tracks
\r
1384 //and return the result
\r
1386 Double_t signDca=-999.;
\r
1387 if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);
\r
1389 //=====Now calculate DCA between both tracks=======
\r
1390 Double_t massE = 0.000511;
\r
1391 Double_t massK = 0.493677;
\r
1393 Double_t vertex[3] = {-999.,-999.,-999}; //vertex
\r
1394 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
\r
1395 GetReader()->GetVertex(vertex); //If only one file, get the vertex from there
\r
1396 //FIXME: Add a check for whether file 2 is PYTHIA or HIJING
\r
1397 //If PYTHIA, then set the vertex from file 2, if not, use the
\r
1398 //vertex from file 1
\r
1399 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);
\r
1402 TVector3 primV(vertex[0],vertex[1],vertex[2]) ;
\r
1404 if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;
\r
1406 AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);
\r
1407 AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);
\r
1409 Double_t bfield[3];
\r
1410 param1->GetBxByBz(bfield);
\r
1411 Double_t bz = param1->GetBz();
\r
1413 Double_t xplane1 = 0.; Double_t xplane2 = 0.;
\r
1414 Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2);
\r
1416 param1->PropagateToBxByBz(xplane1,bfield);
\r
1417 param2->PropagateToBxByBz(xplane2,bfield);
\r
1419 Int_t id1 = 0, id2 = 0;
\r
1420 AliESDv0 bvertex(*param1,id1,*param2,id2);
\r
1421 Double_t vx,vy,vz;
\r
1422 bvertex.GetXYZ(vx,vy,vz);
\r
1426 param1->PxPyPz(emom);
\r
1427 param2->PxPyPz(hmom);
\r
1428 TVector3 emomAtB(emom[0],emom[1],emom[2]);
\r
1429 TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);
\r
1430 TVector3 secvtxpt(vx,vy,vz);
\r
1431 TVector3 decayvector(0,0,0);
\r
1432 decayvector = secvtxpt - primV; //decay vector from PrimVtx
\r
1433 Double_t decaylength = decayvector.Mag();
\r
1435 printf("\t JLK pairDCA = %2.2f\n",pairdca);
\r
1437 if(GetDebug() > 0) {
\r
1438 printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );
\r
1439 printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );
\r
1442 if (masscut<1.1) fhDVMBtagQA1->Fill(pairdca,decaylength);
\r
1444 if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {
\r
1445 TVector3 sumMom = emomAtB+hmomAtB;
\r
1446 Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);
\r
1447 Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);
\r
1448 Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);
\r
1449 Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
\r
1450 Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));
\r
1451 Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();
\r
1453 if (masscut<1.1) fhDVMBtagQA2->Fill(sDca, mass);
\r
1455 if (mass > masscut && massPhot > 0.1) signDca = sDca;
\r
1457 if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);
\r
1458 if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);
\r
1468 //__________________________________________________________________
\r
1469 Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi)
\r
1471 //get signed impact parameter significance of the given AOD track
\r
1472 //for the given jet
\r
1474 Int_t trackIndex = 0;
\r
1475 Int_t ntrk = GetAODCTS()->GetEntriesFast();
\r
1476 for (Int_t k2 =0; k2 < ntrk ; k2++) {
\r
1478 AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);
\r
1479 int id1 = tr->GetID();
\r
1480 int id2 = track2->GetID();
\r
1482 trackIndex = k2;//FIXME: check if GetAODCTS stores tracks in the
\r
1483 //same order of the event
\r
1488 Double_t significance=0;
\r
1489 Double_t maxD = 10000.;
\r
1490 Double_t impPar[] = {0,0};
\r
1491 Double_t ipCov[]={0,0,0};
\r
1492 Double_t ipVec2D[] = {0,0};
\r
1494 AliVEvent* vEvent = (AliVEvent*)GetReader()->GetInputEvent();
\r
1495 if(!vEvent) return -97;
\r
1496 AliVVertex* vv = (AliVVertex*)vEvent->GetPrimaryVertex();
\r
1497 if(!vv) return -98;
\r
1498 AliVTrack* vTrack = (AliVTrack*)vEvent->GetTrack(trackIndex);
\r
1499 if(!vTrack) return -99;
\r
1500 AliESDtrack esdTrack(vTrack);
\r
1501 Double_t bfield[3];
\r
1502 esdTrack.GetBxByBz(bfield);
\r
1503 if(!esdTrack.PropagateToDCABxByBz(vv, bfield, maxD, impPar, ipCov)) return -100;
\r
1504 if(ipCov[0]<0) return -101;
\r
1506 Double_t Pxy[] = {esdTrack.Px(), esdTrack.Py()};
\r
1507 Double_t Txy[] = {esdTrack.Xv(), esdTrack.Yv()};
\r
1508 Double_t Vxy[] = {vv->GetX(), vv->GetY()};
\r
1509 GetImpactParamVect(Pxy, Txy, Vxy, ipVec2D);
\r
1510 Double_t phiIP = TMath::ATan2(ipVec2D[1], ipVec2D[0]) + (TMath::Abs(ipVec2D[1])-ipVec2D[1])/TMath::Abs(ipVec2D[1])*TMath::Pi();
\r
1511 Double_t cosTheta = TMath::Cos(jetPhi - phiIP);
\r
1512 Double_t sign = cosTheta/TMath::Abs(cosTheta);
\r
1513 significance = TMath::Abs(impPar[0])/TMath::Sqrt(ipCov[0])*sign;
\r
1514 printf("\t JLK significance = %2.2f\n",significance);
\r
1515 //ip = fabs(impPar[0]);
\r
1516 fhIPSigBtagQA2->Fill(significance);
\r
1517 return significance;
\r
1520 //__________________________________________________________________
\r
1521 void AliAnaElectron::GetImpactParamVect(Double_t Pxy[2], Double_t Txy[2], Double_t Vxy[2], Double_t IPxy[2])
\r
1523 //px,py: momentum components at the origin of the track; tx, ty:
\r
1524 //origin (x,y) of track; vx, vy: coordinates of primary vertex
\r
1525 // analytical geometry auxiliary variables
\r
1526 Double_t mr = Pxy[1]/Pxy[0]; //angular coeficient of the straight
\r
1527 //line that lies on top of track
\r
1529 Double_t br = Txy[1] - mr*Txy[0]; //linear coeficient of the straight
\r
1530 //line that lies on top of track
\r
1532 Double_t ms = -1./mr; //angular coeficient of the straight line that
\r
1533 //lies on top of the impact parameter line
\r
1534 // Double_t bs = Vxy[1] - ms*Vxy[0]; //linear coeficient of the straight
\r
1535 //line that lies on top of the
\r
1536 //impact parameter line
\r
1537 Double_t xIntersection = (mr*Txy[0] - ms*Vxy[0] + Vxy[1] - Txy[1])/(mr - ms);
\r
1538 Double_t yIntersection = mr*xIntersection + br;
\r
1539 //if(ceil(10000*yIntersection) - ceil(10000*(ms*xIntersection + bs))
\r
1540 //!= 0 )cout<<yIntersection<<", "<<ms*xIntersection + bs<<endl;
\r
1541 IPxy[0] = xIntersection - Vxy[0];
\r
1542 IPxy[1] = yIntersection - Vxy[1];
\r
1546 //__________________________________________________________________
\r
1547 Bool_t AliAnaElectron::PhotonicPrim(const AliAODPWG4Particle* part)
\r
1549 //This method checks the opening angle and invariant mass of
\r
1550 //electron pairs within the AliAODPWG4Particle list to see if
\r
1551 //they are likely to be photonic electrons
\r
1553 Bool_t itIS = kFALSE;
\r
1555 Double_t massE = 0.000511;
\r
1556 Double_t massEta = 0.547;
\r
1557 Double_t massRho0 = 0.770;
\r
1558 Double_t massOmega = 0.782;
\r
1559 Double_t massPhi = 1.020;
\r
1561 Int_t pdg1 = part->GetPdg();
\r
1562 Int_t trackId = part->GetTrackLabel(0);
\r
1563 AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);
\r
1565 if(GetDebug() > 0) printf("AliAnaElectron::PhotonicPrim - can't get the AOD Track from the particle! Skipping the photonic check");
\r
1566 return kFALSE; //Don't proceed because we can't get the track
\r
1569 AliExternalTrackParam *param1 = new AliExternalTrackParam(track);
\r
1571 //Loop on stored AOD electrons and compute the angle differences and Minv
\r
1572 for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {
\r
1573 AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);
\r
1574 Int_t track2Id = part2->GetTrackLabel(0);
\r
1575 if(trackId == track2Id) continue;
\r
1576 Int_t pdg2 = part2->GetPdg();
\r
1577 if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;
\r
1578 if(part2->GetDetector() != fCalorimeter) continue;
\r
1580 //JLK: Check opp. sign pairs only
\r
1581 if(pdg1*pdg2 > 0) continue; //skip same-sign pairs
\r
1583 //propagate to common vertex and check opening angle
\r
1584 AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(track2Id);
\r
1586 if(GetDebug() >0) printf("AliAnaElectron::PhotonicPrim - problem getting the partner track. Continuing on to the next one");
\r
1589 AliExternalTrackParam *param2 = new AliExternalTrackParam(track2);
\r
1590 Int_t id1 = 0, id2 = 0;
\r
1591 AliESDv0 photonVtx(*param1,id1,*param2,id2);
\r
1592 Double_t vx,vy,vz;
\r
1593 photonVtx.GetXYZ(vx,vy,vz);
\r
1595 Double_t p1mom[3];
\r
1596 Double_t p2mom[3];
\r
1597 param1->PxPyPz(p1mom);
\r
1598 param2->PxPyPz(p2mom);
\r
1600 TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);
\r
1601 TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);
\r
1602 TVector3 sumMom = p1momAtB+p2momAtB;
\r
1604 Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);
\r
1605 Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);
\r
1606 Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
\r
1608 Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);
\r
1609 fh1OpeningAngle->Fill(dphi);
\r
1610 fh1MinvPhoton->Fill(mass);
\r
1613 (mass > massEta-0.05 || mass < massEta+0.05) ||
\r
1614 (mass > massRho0-0.05 || mass < massRho0+0.05) ||
\r
1615 (mass > massOmega-0.05 || mass < massOmega+0.05) ||
\r
1616 (mass > massPhi-0.05 || mass < massPhi+0.05))
\r
1619 if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");
\r
1633 //__________________________________________________________________
\r
1634 Bool_t AliAnaElectron::PhotonicV0(Int_t id)
\r
1636 //This method checks to see whether a track that has been flagged as
\r
1637 //an electron was determined to match to a V0 candidate with
\r
1638 //invariant mass consistent with photon conversion
\r
1640 Bool_t itIS = kFALSE;
\r
1642 Double_t massEta = 0.547;
\r
1643 Double_t massRho0 = 0.770;
\r
1644 Double_t massOmega = 0.782;
\r
1645 Double_t massPhi = 1.020;
\r
1648 AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();
\r
1649 int nv0s = aod->GetNumberOfV0s();
\r
1650 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
\r
1651 AliAODv0 *v0 = aod->GetV0(iV0);
\r
1652 if (!v0) continue;
\r
1653 double radius = v0->RadiusV0();
\r
1654 double mass = v0->InvMass2Prongs(0,1,11,11);
\r
1655 if(GetDebug() > 0) {
\r
1656 printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );
\r
1657 printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);
\r
1658 printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );
\r
1660 if (mass < 0.100 ||
\r
1661 (mass > massEta-0.05 || mass < massEta+0.05) ||
\r
1662 (mass > massRho0-0.05 || mass < massRho0+0.05) ||
\r
1663 (mass > massOmega-0.05 || mass < massOmega+0.05) ||
\r
1664 (mass > massPhi-0.05 || mass < massPhi+0.05)) {
\r
1665 if ( id == v0->GetNegID() || id == v0->GetPosID()) {
\r
1667 if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );
\r
1674 //__________________________________________________________________
\r
1675 Bool_t AliAnaElectron::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3])
\r
1677 //Use the Event vertex and AOD track information to get
\r
1678 //a real impact parameter for the track
\r
1679 //Once alice-off gets its act together and fixes the AOD, this
\r
1680 //should become obsolete.
\r
1682 Double_t maxD = 100000.; //max transverse IP
\r
1683 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
\r
1684 AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();
\r
1685 AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();
\r
1686 AliESDtrack esdTrack(track);
\r
1687 Double_t bfield[3];
\r
1688 esdTrack.GetBxByBz(bfield);
\r
1689 Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov);
\r
1690 printf("\t JLK impPar = %2.2f\n",impPar[0]);
\r
1698 //__________________________________________________________________
\r
1699 Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type)
\r
1701 //Check this track to see if it is also tagged as an electron in the
\r
1702 //AliAODPWG4Particle list and if it is non-photonic
\r
1704 Bool_t pass = kFALSE;
\r
1706 Int_t trackId = track->GetID(); //get the index in the reader
\r
1708 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
\r
1709 if(GetDebug() > 3) printf("AliAnaElectron::CheckTrack() - aod branch entries %d\n", naod);
\r
1710 for(Int_t iaod = 0; iaod < naod ; iaod++){
\r
1711 AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
\r
1712 Int_t label = ele->GetTrackLabel(0);
\r
1713 if(label != trackId) continue; //skip to the next one if they don't match
\r
1715 if(strcmp(type,"DVM")==0) {
\r
1716 if(ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag1) ||
\r
1717 ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag2))
\r
1720 } else if (strcmp(type,"NPE")==0) {
\r
1722 Bool_t photonic = kFALSE;
\r
1723 Bool_t photonic1 = kFALSE;
\r
1724 photonic1 = PhotonicPrim(ele); //check against primaries
\r
1725 Bool_t photonic2 = kFALSE;
\r
1726 photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s
\r
1727 if(photonic1 || photonic2) photonic = kTRUE;
\r
1729 if(!photonic) pass = kTRUE;
\r
1740 //__________________________________________________________________
\r
1741 Double_t AliAnaElectron::GetBParentPt(Int_t ipart)
\r
1743 //return MC B parent pt
\r
1744 if(GetReader()->ReadStack()) { //only done if we have the stack
\r
1745 AliStack* stack = GetMCStack();
\r
1747 printf("Problem getting stack\n");
\r
1750 TParticle* prim = stack->Particle(ipart);
\r
1751 if(prim->GetMother(0)>=0) {
\r
1753 TParticle *parent = stack->Particle(prim->GetMother(0));
\r
1754 if(parent) mpdg = parent->GetPdgCode();
\r
1756 if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) ||
\r
1757 (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )
\r
1758 return parent->Pt();
\r
1760 } else if(GetReader()->ReadAODMCParticles()){
\r
1761 AliAODMCParticle* prim = (AliAODMCParticle*)GetMCParticle(ipart);
\r
1762 if(prim->GetMother()>=0) {
\r
1764 AliAODMCParticle* parent = (AliAODMCParticle*)GetMCParticle(prim->GetMother());
\r
1765 if(parent) mpdg = parent->GetPdgCode();
\r
1766 if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) ||
\r
1767 (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )
\r
1768 return parent->Pt();
\r
1774 //__________________________________________________________________
\r
1775 Int_t AliAnaElectron::GetMCSource(Int_t tag)
\r
1777 //For determining how to classify electrons using MC info
\r
1778 //the number returned is the bin along one axis of 2-d histograms in
\r
1779 //which to fill this electron
\r
1782 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;
\r
1784 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) {
\r
1786 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 1;
\r
1788 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC)
\r
1789 && !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 2;
\r
1790 //Charm from bottom
\r
1791 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB)) return 3;
\r
1793 //else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;
\r
1795 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)
\r
1796 || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay)
\r
1797 || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) return 5;
\r
1799 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay)
\r
1800 || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay)) return 6;
\r
1805 //Misidentified electron
\r
1811 //__________________________________________________________________
\r
1812 Int_t AliAnaElectron::GetNumAODMCParticles()
\r
1814 //Get the number of AliAODMCParticles, if any
\r
1817 TClonesArray * mcparticles0 = 0x0;
\r
1818 TClonesArray * mcparticles1 = 0x0;
\r
1820 if(GetReader()->ReadAODMCParticles()){
\r
1821 //Get the list of MC particles
\r
1823 mcparticles0 = GetReader()->GetAODMCParticles(0);
\r
1824 if(!mcparticles0 && GetDebug() > 0) {
\r
1825 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
\r
1827 // if(GetReader()->GetSecondInputAODTree()){
\r
1828 // mcparticles1 = GetReader()->GetAODMCParticles(1);
\r
1829 // if(!mcparticles1 && GetDebug() > 0) {
\r
1830 // printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
\r
1834 Int_t npart0 = mcparticles0->GetEntriesFast();
\r
1836 if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();
\r
1837 Int_t npart = npart0+npart1;
\r
1844 //__________________________________________________________________
\r
1845 AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t ipart)
\r
1847 //Get the MC particle at position ipart
\r
1849 AliAODMCParticle* aodprimary = 0x0;
\r
1850 TClonesArray * mcparticles0 = 0x0;
\r
1851 TClonesArray * mcparticles1 = 0x0;
\r
1853 if(GetReader()->ReadAODMCParticles()){
\r
1854 //Get the list of MC particles
\r
1855 mcparticles0 = GetReader()->GetAODMCParticles(0);
\r
1856 if(!mcparticles0 && GetDebug() > 0) {
\r
1857 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
\r
1859 // if(GetReader()->GetSecondInputAODTree()){
\r
1860 // mcparticles1 = GetReader()->GetAODMCParticles(1);
\r
1861 // if(!mcparticles1 && GetDebug() > 0) {
\r
1862 // printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
\r
1866 Int_t npart0 = mcparticles0->GetEntriesFast();
\r
1868 if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();
\r
1869 if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);
\r
1870 else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);
\r
1872 printf("AliAnaElectron::GetMCParticle() *** no primary ***: label %d \n", ipart);
\r
1877 printf("AliAnaElectron::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");
\r
1879 return aodprimary;
\r
1883 //__________________________________________________________________
\r
1884 Bool_t AliAnaElectron::IsMcBJet(Double_t jeta, Double_t jphi)
\r
1886 //Check the jet eta,phi against that of the b-quark
\r
1887 //to decide whether it is an MC B-jet
\r
1888 Bool_t bjet=kFALSE;
\r
1890 // printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() );
\r
1892 AliStack* stack = 0x0;
\r
1894 for(Int_t ipart = 0; ipart < 100; ipart++) {
\r
1896 Double_t pphi = -999.;
\r
1897 Double_t peta = -999.;
\r
1899 if(GetReader()->ReadStack()) {
\r
1900 stack = GetMCStack();
\r
1902 printf("AliAnaElectron::IsMCBJet() *** no stack ***: \n");
\r
1905 TParticle* primary = stack->Particle(ipart);
\r
1906 if (!primary) continue;
\r
1907 pdg = primary->GetPdgCode();
\r
1908 pphi = primary->Phi();
\r
1909 peta = primary->Eta();
\r
1910 } else if(GetReader()->ReadAODMCParticles()) {
\r
1911 AliAODMCParticle* aodprimary = GetMCParticle(ipart);
\r
1912 if(!aodprimary) continue;
\r
1913 pdg = aodprimary->GetPdgCode();
\r
1914 pphi = aodprimary->Phi();
\r
1915 peta = aodprimary->Eta();
\r
1917 if ( TMath::Abs(pdg) != 5) continue;
\r
1919 // printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());
\r
1920 Double_t dphi = jphi - pphi;
\r
1921 Double_t deta = jeta - peta;
\r
1922 Double_t dr = sqrt(deta*deta + dphi*dphi);
\r
1926 //printf("MTH: **** found matching MC-Bjet: PDG=%d, pt=%f,dr=%f \n", pdgcode, primary->Pt(),dr );
\r
1934 //__________________________________________________________________
\r
1935 Bool_t AliAnaElectron::IsMcDJet(Double_t jeta, Double_t jphi)
\r
1937 //Check if this jet is a charm jet
\r
1938 Bool_t cjet=kFALSE;
\r
1940 AliStack* stack = 0x0;
\r
1942 for(Int_t ipart = 0; ipart < 100; ipart++) {
\r
1944 Double_t pphi = -999.;
\r
1945 Double_t peta = -999.;
\r
1947 if(GetReader()->ReadStack()) {
\r
1948 stack = GetMCStack();
\r
1950 printf("AliAnaElectron::IsMCDJet() *** no stack ***: \n");
\r
1953 TParticle* primary = stack->Particle(ipart);
\r
1954 if (!primary) continue;
\r
1955 pdg = primary->GetPdgCode();
\r
1956 pphi = primary->Phi();
\r
1957 peta = primary->Eta();
\r
1958 } else if(GetReader()->ReadAODMCParticles()) {
\r
1959 AliAODMCParticle* aodprimary = GetMCParticle(ipart);
\r
1960 if(!aodprimary) continue;
\r
1961 pdg = aodprimary->GetPdgCode();
\r
1962 pphi = aodprimary->Phi();
\r
1963 peta = aodprimary->Eta();
\r
1966 if ( TMath::Abs(pdg) != 4) continue;
\r
1968 Double_t dphi = jphi - pphi;
\r
1969 Double_t deta = jeta - peta;
\r
1970 Double_t dr = sqrt(deta*deta + dphi*dphi);
\r
1982 //__________________________________________________________________
\r
1983 void AliAnaElectron::Print(const Option_t * opt) const
\r
1985 //Print some relevant parameters set for the analysis
\r
1990 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
\r
1991 AliAnaPartCorrBaseClass::Print(" ");
\r
1993 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
\r
1994 printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax);
\r
1995 printf("residual cut = %f\n",fResidualCut);
\r
1996 printf("---DVM Btagging\n");
\r
1997 printf("max IP-cut (e,h) = %f\n",fImpactCut);
\r
1998 printf("min ITS-hits = %d\n",fITSCut);
\r
1999 printf("max dR (e,h) = %f\n",fDrCut);
\r
2000 printf("max pairDCA = %f\n",fPairDcaCut);
\r
2001 printf("max decaylength = %f\n",fDecayLenCut);
\r
2002 printf("min Associated Pt = %f\n",fAssocPtCut);
\r
2003 printf("---IPSig Btagging\n");
\r
2004 printf("min tag track = %d\n",fNTagTrkCut);
\r
2005 printf("min IP significance = %f\n",fIPSigCut);
\r
2010 //________________________________________________________________________
\r
2011 void AliAnaElectron::ReadHistograms(TList* /* outputList */)
\r
2013 // Needed when Terminate is executed in distributed environment
\r
2014 // Refill analysis histograms of this class with corresponding
\r
2015 // histograms in output list.
\r
2017 // Histograms of this analsys are kept in the same list as other
\r
2018 // analysis, recover the position of
\r
2019 // the first one and then add the next
\r
2020 //Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));
\r
2022 //Read histograms, must be in the same order as in
\r
2023 //GetCreateOutputObject.
\r
2024 //fh1pOverE = (TH1F *) outputList->At(index);
\r
2025 //fh1dR = (TH1F *) outputList->At(index++);
\r
2026 //fh2EledEdx = (TH2F *) outputList->At(index++);
\r
2027 //fh2MatchdEdx = (TH2F *) outputList->At(index++);
\r
2031 //__________________________________________________________________
\r
2032 void AliAnaElectron::Terminate(TList* outputList)
\r
2035 //Do some plots to end
\r
2036 //Recover histograms from output histograms list, needed for
\r
2037 //distributed analysis.
\r
2038 //ReadHistograms(outputList);
\r
2040 printf(" AliAnaElectron::Terminate() *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;
\r
git://git.uio.no / u / mrichter / AliRoot.git / blob