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