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