1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //_________________________________________________________________________
19 // Class for the electron identification.
20 // Clusters from EMCAL matched to tracks
21 // and kept in the AOD. Few histograms produced.
23 // -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)
24 //////////////////////////////////////////////////////////////////////////////
26 // --- ROOT system ---
28 #include <TParticle.h>
30 #include <TClonesArray.h>
31 #include <TObjString.h>
32 //#include <Riostream.h>
34 // --- Analysis system ---
35 #include "AliAnaElectron.h"
36 #include "AliCaloTrackReader.h"
37 #include "AliMCAnalysisUtils.h"
38 #include "AliAODCaloCluster.h"
39 #include "AliFidutialCut.h"
40 #include "AliAODTrack.h"
41 #include "AliAODPid.h"
42 #include "AliCaloPID.h"
43 #include "AliAODMCParticle.h"
45 #include "AliExternalTrackParam.h"
48 ClassImp(AliAnaElectron)
50 //____________________________________________________________________________
51 AliAnaElectron::AliAnaElectron()
52 : AliAnaPartCorrBaseClass(),fCalorimeter(""),
53 fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),
54 fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),
55 fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),
59 fh1pOverE(0),fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),
60 fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),
61 fh2TrackPVsClusterE(0),fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),
62 //Photonic electron checks
63 fh1OpeningAngle(0),fh1MinvPhoton(0),
65 fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
66 fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),
67 fhPtPE(0),fhPhiPE(0),fhEtaPE(0),
68 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
69 fhPtBottom(0),fhPhiBottom(0),fhEtaBottom(0),
70 fhPtCharm(0),fhPhiCharm(0),fhEtaCharm(0),
71 fhPtCFromB(0),fhPhiCFromB(0),fhEtaCFromB(0),
72 fhPtDalitz(0),fhPhiDalitz(0),fhEtaDalitz(0),
73 fhPtWDecay(0),fhPhiWDecay(0),fhEtaWDecay(0),
74 fhPtZDecay(0),fhPhiZDecay(0),fhEtaZDecay(0),
75 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
76 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
78 fhBtagCut1(0),fhBtagCut2(0),fhBtagCut3(0),
84 //Initialize parameters
89 //____________________________________________________________________________
90 AliAnaElectron::AliAnaElectron(const AliAnaElectron & g)
91 : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),
92 fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),fResidualCut(g.fResidualCut),
93 fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),
94 fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),
95 fWriteNtuple(g.fWriteNtuple),
97 fEleNtuple(g.fEleNtuple),
98 fh1pOverE(g.fh1pOverE),fh1dR(g.fh1dR),
99 fh2EledEdx(g.fh2EledEdx),fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),
100 fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),
101 fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),
102 fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),
103 //Photonic electron checks
104 fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),
106 fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron),
107 fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE),
108 fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE),
109 fhPtConversion(g.fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
110 fhPtBottom(g.fhPtBottom),fhPhiBottom(g.fhPhiBottom),fhEtaBottom(g.fhEtaBottom),
111 fhPtCharm(g.fhPtCharm),fhPhiCharm(g.fhPhiCharm),fhEtaCharm(g.fhEtaCharm),
112 fhPtCFromB(g.fhPtCFromB),fhPhiCFromB(g.fhPhiCFromB),fhEtaCFromB(g.fhEtaCFromB),
113 fhPtDalitz(g.fhPtDalitz),fhPhiDalitz(g.fhPhiDalitz),fhEtaDalitz(g.fhEtaDalitz),
114 fhPtWDecay(g.fhPtWDecay),fhPhiWDecay(g.fhPhiWDecay),fhEtaWDecay(g.fhEtaWDecay),
115 fhPtZDecay(g.fhPtZDecay),fhPhiZDecay(g.fhPhiZDecay),fhEtaZDecay(g.fhEtaZDecay),
116 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
117 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
119 fhBtagCut1(g.fhBtagCut1),fhBtagCut2(g.fhBtagCut2),fhBtagCut3(g.fhBtagCut3),
121 fMCEleNtuple(g.fMCEleNtuple)
127 //_________________________________________________________________________
128 AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)
130 // assignment operator
132 if(&g == this) return *this;
133 fCalorimeter = g.fCalorimeter;
134 fpOverEmin = g.fpOverEmin;
135 fpOverEmax = g.fpOverEmax;
136 fResidualCut = g.fResidualCut;
138 fPairDcaCut = g.fPairDcaCut;
139 fDecayLenCut = g.fDecayLenCut;
140 fImpactCut = g.fImpactCut;
141 fAssocPtCut = g.fAssocPtCut;
142 fMassCut = g.fMassCut;
143 fSdcaCut = g.fSdcaCut;
145 fWriteNtuple = g.fWriteNtuple;
146 fEleNtuple = g.fEleNtuple;
147 fh1pOverE = g.fh1pOverE;
149 fh2EledEdx = g.fh2EledEdx;
150 fh2MatchdEdx = g.fh2MatchdEdx;
151 fh2dEtadPhi = g.fh2dEtadPhi;
152 fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;
153 fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;
154 fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;
155 fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;
156 fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;
157 fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;
158 fh1OpeningAngle = g.fh1OpeningAngle;
159 fh1MinvPhoton = g.fh1MinvPhoton;
160 fhPtElectron = g.fhPtElectron;
161 fhPhiElectron = g.fhPhiElectron;
162 fhEtaElectron = g.fhEtaElectron;
164 fhPhiNPE = g.fhPhiNPE;
165 fhEtaNPE = g.fhEtaNPE;
169 fhPtConversion = g.fhPtConversion;
170 fhPhiConversion = g.fhPhiConversion;
171 fhEtaConversion = g.fhEtaConversion;
172 fhPtBottom = g.fhPtBottom;
173 fhPhiBottom = g.fhPhiBottom;
174 fhEtaBottom = g.fhEtaBottom;
175 fhPtCharm = g.fhPtCharm;
176 fhPhiCharm = g.fhPhiCharm;
177 fhEtaCharm = g.fhEtaCharm;
178 fhPtCFromB = g.fhPtCFromB;
179 fhPhiCFromB = g.fhPhiCFromB;
180 fhEtaCFromB = g.fhEtaCFromB;
181 fhPtDalitz = g.fhPtDalitz;
182 fhPhiDalitz = g.fhPhiDalitz;
183 fhEtaDalitz = g.fhEtaDalitz;
184 fhPtWDecay = g.fhPtWDecay;
185 fhPhiWDecay = g.fhPhiWDecay;
186 fhEtaWDecay = g.fhEtaWDecay;
187 fhPtZDecay = g.fhPtZDecay;
188 fhPhiZDecay = g.fhPhiZDecay;
189 fhEtaZDecay = g.fhEtaZDecay;
190 fhPtPrompt = g.fhPtPrompt;
191 fhPhiPrompt = g.fhPhiPrompt;
192 fhEtaPrompt = g.fhEtaPrompt;
193 fhPtUnknown = g.fhPtUnknown;
194 fhPhiUnknown = g.fhPhiUnknown;
195 fhEtaUnknown = g.fhEtaUnknown;
196 fMCEleNtuple = g.fMCEleNtuple;
199 fhBtagCut1 = g.fhBtagCut1;
200 fhBtagCut2 = g.fhBtagCut2;
201 fhBtagCut3 = g.fhBtagCut3;
207 //____________________________________________________________________________
208 AliAnaElectron::~AliAnaElectron()
215 //________________________________________________________________________
216 TList * AliAnaElectron::GetCreateOutputObjects()
218 // Create histograms to be saved in output file and
219 // store them in outputContainer
220 TList * outputContainer = new TList() ;
221 outputContainer->SetName("ElectronHistos") ;
223 //created ele ntuple for further analysis
225 fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");
226 outputContainer->Add(fEleNtuple) ;
229 Int_t nptbins = GetHistoNPtBins();
230 Int_t nphibins = GetHistoNPhiBins();
231 Int_t netabins = GetHistoNEtaBins();
232 Float_t ptmax = GetHistoPtMax();
233 Float_t phimax = GetHistoPhiMax();
234 Float_t etamax = GetHistoEtaMax();
235 Float_t ptmin = GetHistoPtMin();
236 Float_t phimin = GetHistoPhiMin();
237 Float_t etamin = GetHistoEtaMin();
239 fh1pOverE = new TH1F("h1pOverE","EMCAL-TRACK matches p/E",100,0.,10.);
240 fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi());
241 fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.);
242 fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.);
243 fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
244 fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
245 fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
247 fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
248 fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
249 fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax);
250 fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax);
252 outputContainer->Add(fh1pOverE) ;
253 outputContainer->Add(fh1dR) ;
254 outputContainer->Add(fh2EledEdx) ;
255 outputContainer->Add(fh2MatchdEdx) ;
256 outputContainer->Add(fh2dEtadPhi) ;
257 outputContainer->Add(fh2dEtadPhiMatched) ;
258 outputContainer->Add(fh2dEtadPhiUnmatched) ;
259 outputContainer->Add(fh2TrackPVsClusterE) ;
260 outputContainer->Add(fh2TrackPtVsClusterE) ;
261 outputContainer->Add(fh2TrackPhiVsClusterPhi) ;
262 outputContainer->Add(fh2TrackEtaVsClusterEta) ;
264 //photonic electron checks
265 fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between electron pairs",100,0.,TMath::Pi());
266 fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of electron pairs",100,0.,2.);
268 outputContainer->Add(fh1OpeningAngle);
269 outputContainer->Add(fh1MinvPhoton);
271 fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);
272 fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
273 fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
274 fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);
275 fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
276 fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
277 fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);
278 fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
279 fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
281 outputContainer->Add(fhPtElectron) ;
282 outputContainer->Add(fhPhiElectron) ;
283 outputContainer->Add(fhEtaElectron) ;
284 outputContainer->Add(fhPtNPE) ;
285 outputContainer->Add(fhPhiNPE) ;
286 outputContainer->Add(fhEtaNPE) ;
287 outputContainer->Add(fhPtPE) ;
288 outputContainer->Add(fhPhiPE) ;
289 outputContainer->Add(fhEtaPE) ;
292 fhBtagCut1 = new TH2F("hbtag_cut1","B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);
293 fhBtagCut2 = new TH2F("hbtag_cut2","B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);
294 fhBtagCut3 = new TH2F("hbtag_cut3","B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);
296 outputContainer->Add(fhBtagCut1) ;
297 outputContainer->Add(fhBtagCut2) ;
298 outputContainer->Add(fhBtagCut3) ;
302 fhPtConversion = new TH1F("hPtConversion","Conversion electron pT",nptbins,ptmin,ptmax);
303 fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
304 fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
305 fhPtBottom = new TH1F("hPtBottom","Bottom electron pT",nptbins,ptmin,ptmax);
306 fhPhiBottom = new TH2F("hPhiBottom","Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
307 fhEtaBottom = new TH2F("hEtaBottom","Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
308 fhPtCharm = new TH1F("hPtCharm","Charm electron pT",nptbins,ptmin,ptmax);
309 fhPhiCharm = new TH2F("hPhiCharm","Charm Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
310 fhEtaCharm = new TH2F("hEtaCharm","Charm Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
311 fhPtCFromB = new TH1F("hPtCFromB","Charm from Bottom electron pT",nptbins,ptmin,ptmax);
312 fhPhiCFromB = new TH2F("hPhiCFromB","Charm from Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
313 fhEtaCFromB = new TH2F("hEtaCFromB","Charm from Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
314 fhPtDalitz = new TH1F("hPtDalitz","Dalitz electron pT",nptbins,ptmin,ptmax);
315 fhPhiDalitz = new TH2F("hPhiDalitz","Dalitz Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
316 fhEtaDalitz = new TH2F("hEtaDalitz","Dalitz Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
317 fhPtWDecay = new TH1F("hPtWDecay","W-boson Electron pT",nptbins,ptmin,ptmax);
318 fhPhiWDecay = new TH2F("hPhiWDecay","W-boson electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
319 fhEtaWDecay = new TH2F("hEtaWDecay","W-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
320 fhPtZDecay = new TH1F("hPtZDecay","Z-boson electron pT",nptbins,ptmin,ptmax);
321 fhPhiZDecay = new TH2F("hPhiZDecay","Z-boson Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
322 fhEtaZDecay = new TH2F("hEtaZDecay","Z-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
323 fhPtPrompt = new TH1F("hPtPrompt","Prompt electron pT",nptbins,ptmin,ptmax);
324 fhPhiPrompt = new TH2F("hPhiPrompt","Prompt Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
325 fhEtaPrompt = new TH2F("hEtaPrompt","Prompt Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
326 fhPtUnknown = new TH1F("hPtUnknown","Unknown electron pT",nptbins,ptmin,ptmax);
327 fhPhiUnknown = new TH2F("hPhiUnknown","Unknown Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
328 fhEtaUnknown = new TH2F("hEtaUnknown","Unknown Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
330 outputContainer->Add(fhPtConversion);
331 outputContainer->Add(fhPhiConversion);
332 outputContainer->Add(fhEtaConversion);
333 outputContainer->Add(fhPtBottom);
334 outputContainer->Add(fhPhiBottom);
335 outputContainer->Add(fhEtaBottom);
336 outputContainer->Add(fhPtCharm);
337 outputContainer->Add(fhPhiCharm);
338 outputContainer->Add(fhEtaCharm);
339 outputContainer->Add(fhPtCFromB);
340 outputContainer->Add(fhPhiCFromB);
341 outputContainer->Add(fhEtaCFromB);
342 outputContainer->Add(fhPtDalitz);
343 outputContainer->Add(fhPhiDalitz);
344 outputContainer->Add(fhEtaDalitz);
345 outputContainer->Add(fhPtWDecay);
346 outputContainer->Add(fhPhiWDecay);
347 outputContainer->Add(fhEtaWDecay);
348 outputContainer->Add(fhPtZDecay);
349 outputContainer->Add(fhPhiZDecay);
350 outputContainer->Add(fhEtaZDecay);
351 outputContainer->Add(fhPtPrompt);
352 outputContainer->Add(fhPhiPrompt);
353 outputContainer->Add(fhEtaPrompt);
354 outputContainer->Add(fhPtUnknown);
355 outputContainer->Add(fhPhiUnknown);
356 outputContainer->Add(fhEtaUnknown);
358 //created ele ntuple for further analysis
360 fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");
361 outputContainer->Add(fMCEleNtuple) ;
366 //Save parameters used for analysis
367 TString parList ; //this will be list of parameters used for this analysis.
370 sprintf(onePar,"--- AliAnaElectron ---\n") ;
372 sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ;
374 sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ;
376 sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ;
378 sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ;
380 sprintf(onePar,"---Btagging\n");
382 sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut);
384 sprintf(onePar,"min ITS-hits: %d\n",fITSCut);
386 sprintf(onePar,"max dR (e,h): %f\n",fDrCut);
388 sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut);
390 sprintf(onePar,"max decaylength: %f\n",fDecayLenCut);
392 sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut);
395 //Get parameters set in base class.
396 parList += GetBaseParametersList() ;
398 //Get parameters set in FidutialCut class (not available yet)
399 //parlist += GetFidCut()->GetFidCutParametersList()
401 TObjString *oString= new TObjString(parList) ;
402 outputContainer->Add(oString);
404 return outputContainer ;
408 //____________________________________________________________________________
409 void AliAnaElectron::Init()
412 //do some initialization
413 if(fCalorimeter == "PHOS") {
414 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");
415 fCalorimeter = "EMCAL";
417 if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
418 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");
425 //____________________________________________________________________________
426 void AliAnaElectron::InitParameters()
429 //Initialize the parameters of the analysis.
430 SetOutputAODClassName("AliAODPWG4Particle");
431 SetOutputAODName("PWG4Particle");
433 AddToHistogramsName("AnaElectron_");
435 fCalorimeter = "EMCAL" ;
451 //__________________________________________________________________
452 void AliAnaElectron::MakeAnalysisFillAOD()
455 // Do analysis and fill aods with electron candidates
456 // These AODs will be used to do subsequent histogram filling
458 // Also fill some QA histograms
461 TObjArray *cl = new TObjArray();
463 Double_t bfield = 0.;
464 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
466 //Select the calorimeter of the electron
467 if(fCalorimeter != "EMCAL") {
468 printf("This class not yet implemented for PHOS\n");
473 ////////////////////////////////////////////////
474 //Start from tracks and get associated clusters
475 ////////////////////////////////////////////////
476 if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
477 Int_t ntracks = GetAODCTS()->GetEntriesFast();
479 printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
481 //Unfortunately, AliAODTracks don't have associated EMCAL clusters.
482 //we have to redo track-matching, I guess
483 Int_t iCluster = -999;
484 Int_t bt = 0; //counter for event b-tags
486 for (Int_t itrk = 0; itrk < ntracks; itrk++) {////////////// track loop
487 iCluster = -999; //start with no match
488 AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;
489 AliAODPid* pid = (AliAODPid*) track->GetDetPid();
492 pid->GetEMCALPosition(emcpos);
494 pid->GetEMCALMomentum(emcmom);
496 TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);
497 TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);
498 Double_t tphi = pos.Phi();
499 Double_t teta = pos.Eta();
500 Double_t tmom = mom.Mag();
502 TLorentzVector mom2(mom,0.);
503 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter) ;
504 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);
505 if(mom.Pt() > GetMinPt() && in) {
507 Double_t dEdx = pid->GetTPCsignal();
509 //NOTE: As of 02-Sep-2009, the XYZAtDCA methods of AOD do not
510 //work, but it is possible to get the position of a track at
511 //closest approach to the vertex from the GetPosition method
513 //track->XYZAtDCA(xyz);
514 Bool_t isNotDCA = track->GetPosition(xyz);
515 if(isNotDCA) printf("##Problem getting impact parameter!\n");
516 //printf("\tTRACK POSITION AT DCA: %2.2f,%2.2f,%2.2f\n",xyz[0],xyz[1],xyz[2]);
517 Double_t xy = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
520 Int_t ntot = cl->GetEntriesFast();
522 Double_t pOverE = -999.;
524 Bool_t isElectron = kFALSE;
525 //For tracks in EMCAL acceptance, pair them with all clusters
526 //and fill the dEta vs dPhi for these pairs:
527 for(Int_t iclus = 0; iclus < ntot; iclus++) {
528 AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));
532 clus->GetPosition(x);
533 TVector3 cluspos(x[0],x[1],x[2]);
534 Double_t deta = teta - cluspos.Eta();
535 Double_t dphi = tphi - cluspos.Phi();
536 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
537 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
538 fh2dEtadPhi->Fill(deta,dphi);
539 fh2TrackPVsClusterE->Fill(clus->E(),track->P());
540 fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());
541 fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());
542 fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());
544 res = sqrt(dphi*dphi + deta*deta);
546 fh2dEtadPhiMatched->Fill(deta,dphi);
548 /////////////////////////////////
549 //Perform electron cut analysis//
550 /////////////////////////////////
552 if(res < fResidualCut) {
559 //Input from second AOD?
561 if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;
562 tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);
564 //Do you want the cluster or the track label?
566 if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;
567 cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(0),GetReader(),input);
571 fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,track->GetMostProbablePID(),xy,z);
574 fh2MatchdEdx->Fill(track->P(),dEdx);
576 Double_t energy = clus->E();
577 if(energy > 0) pOverE = tmom/energy;
578 fh1pOverE->Fill(pOverE);
580 Int_t mult = clus->GetNCells();
581 if(mult < 2 && GetDebug() > 0) printf("Single digit cluster.\n");
583 //////////////////////////////
584 //Electron cuts happen here!//
585 //////////////////////////////
586 if(pOverE > fpOverEmin && pOverE < fpOverEmax) isElectron = kTRUE;
588 fh2dEtadPhiUnmatched->Fill(deta,dphi);
593 ///////////////////////////
594 //Fill AOD with electrons//
595 ///////////////////////////
599 if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");
600 Int_t btag = GetBtag(track); bt += btag;
602 fh2EledEdx->Fill(track->P(),dEdx);
604 Double_t eMass = 0.511/1000; //mass in GeV
605 Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);
606 AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);
607 tr.SetLabel(track->GetLabel());
608 tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters
609 tr.SetTrackLabel(track->GetID(),-1); //sets the indices of the original tracks
610 tr.SetDetector(fCalorimeter);
611 if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);
612 //Make this preserve sign of particle
613 if(track->Charge() < 0) tr.SetPdg(11); //electron is 11
614 else tr.SetPdg(-11); //positron is -11
617 //Play with the MC stack if available
618 //Check origin of the candidates
621 //FIXME: Need to re-think this for track-oriented analysis
622 //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?
623 tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));
625 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());
626 }//Work with stack also
630 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());
632 }//pt, fiducial selection
635 //FIXME: Should we also check from the calocluster side, just in
638 if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");
639 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs \n");
643 //__________________________________________________________________
644 void AliAnaElectron::MakeAnalysisFillHistograms()
646 //Do analysis and fill histograms
648 AliStack * stack = 0x0;
649 TParticle * primary = 0x0;
650 TClonesArray * mcparticles0 = 0x0;
651 TClonesArray * mcparticles1 = 0x0;
652 AliAODMCParticle * aodprimary = 0x0;
655 if(GetReader()->ReadStack()){
656 stack = GetMCStack() ;
659 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");
662 else if(GetReader()->ReadAODMCParticles()){
663 //Get the list of MC particles
664 mcparticles0 = GetReader()->GetAODMCParticles(0);
665 if(!mcparticles0 && GetDebug() > 0) {
666 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
668 if(GetReader()->GetSecondInputAODTree()){
669 mcparticles1 = GetReader()->GetAODMCParticles(1);
670 if(!mcparticles1 && GetDebug() > 0) {
671 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
678 //Loop on stored AOD electrons
679 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
680 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
682 for(Int_t iaod = 0; iaod < naod ; iaod++){
683 AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
684 Int_t pdg = ele->GetPdg();
687 printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;
689 if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue;
690 if(ele->GetDetector() != fCalorimeter) continue;
693 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;
696 //Filter for photonic electrons based on opening angle and Minv
697 //cuts, also fill histograms
698 Bool_t photonic = kFALSE;
699 photonic = IsItPhotonic(ele);
701 //Fill electron histograms
702 Float_t ptele = ele->Pt();
703 Float_t phiele = ele->Phi();
704 Float_t etaele = ele->Eta();
706 fhPtElectron ->Fill(ptele);
707 fhPhiElectron ->Fill(ptele,phiele);
708 fhEtaElectron ->Fill(ptele,etaele);
712 fhPhiPE->Fill(ptele,phiele);
713 fhEtaPE->Fill(ptele,etaele);
715 fhPtNPE->Fill(ptele);
716 fhPhiNPE->Fill(ptele,phiele);
717 fhEtaNPE->Fill(ptele,etaele);
721 Int_t tag = ele->GetTag();
722 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)){
723 fhPtConversion ->Fill(ptele);
724 fhPhiConversion ->Fill(ptele,phiele);
725 fhEtaConversion ->Fill(ptele,etaele);
727 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB))
729 printf("\t\tTAG VALUE = %d\n",tag);
730 fhPtBottom ->Fill(ptele);
731 fhPhiBottom ->Fill(ptele,phiele);
732 fhEtaBottom ->Fill(ptele,etaele);
734 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC))
736 fhPtCharm ->Fill(ptele);
737 fhPhiCharm ->Fill(ptele,phiele);
738 fhEtaCharm ->Fill(ptele,etaele);
740 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB))
742 fhPtCFromB ->Fill(ptele);
743 fhPhiCFromB ->Fill(ptele,phiele);
744 fhEtaCFromB ->Fill(ptele,etaele);
746 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
748 fhPtDalitz ->Fill(ptele);
749 fhPhiDalitz ->Fill(ptele,phiele);
750 fhEtaDalitz ->Fill(ptele,etaele);
752 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay))
754 fhPtWDecay ->Fill(ptele);
755 fhPhiWDecay ->Fill(ptele,phiele);
756 fhEtaWDecay ->Fill(ptele,etaele);
758 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay))
760 fhPtZDecay ->Fill(ptele);
761 fhPhiZDecay ->Fill(ptele,phiele);
762 fhEtaZDecay ->Fill(ptele,etaele);
764 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
766 fhPtPrompt ->Fill(ptele);
767 fhPhiPrompt ->Fill(ptele,phiele);
768 fhEtaPrompt ->Fill(ptele,etaele);
771 fhPtUnknown ->Fill(ptele);
772 fhPhiUnknown ->Fill(ptele,phiele);
773 fhEtaUnknown ->Fill(ptele,etaele);
775 }//Histograms with MC
779 ////////////////////////////////////////////////////////
780 //Fill histograms of pure MC kinematics from the stack//
781 ////////////////////////////////////////////////////////
783 if(GetReader()->ReadStack()) {
784 for(Int_t ipart = 0; ipart < stack->GetNtrack(); ipart++) {
785 primary = stack->Particle(ipart);
786 Int_t pdgcode = primary->GetPdgCode();
787 //we only care about electrons
788 if(TMath::Abs(pdgcode) != 11) continue;
789 //we only want TRACKABLE electrons (TPC 85-250cm)
790 if(primary->R() > 200.) continue;
791 //Ignore low pt electrons
792 if(primary->Pt() < 0.2) continue;
794 //find out what the ancestry of this electron is
797 mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);
801 fMCEleNtuple->Fill(mctag,primary->Pt(),primary->Phi(),primary->Eta(),primary->Vx(),primary->Vy(),primary->Vz());
806 } else if(GetReader()->ReadAODMCParticles()) {
807 Int_t npart0 = mcparticles0->GetEntriesFast();
809 if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();
810 Int_t npart = npart0+npart1;
811 for(Int_t ipart = 0; ipart < npart; ipart++) {
812 if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);
813 else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);
815 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", ipart);
818 Int_t pdgcode = aodprimary->GetPdgCode();
819 //we only care about electrons
820 if(TMath::Abs(pdgcode) != 11) continue;
821 //we only want TRACKABLE electrons (TPC 85-250cm)
822 Double_t radius = TMath::Sqrt(aodprimary->Xv()*aodprimary->Xv() + aodprimary->Yv()*aodprimary->Yv());
823 if(radius > 200.) continue;
825 if(aodprimary->Pt() < 0.2) continue;
827 //find out what the ancestry of this electron is
831 if(ipart > npart0) { ival -= npart0; input = 1;}
832 mctag = GetMCAnalysisUtils()->CheckOrigin(ival,GetReader(),input);
836 fMCEleNtuple->Fill(mctag,aodprimary->Pt(),aodprimary->Phi(),aodprimary->Eta(),aodprimary->Xv(),aodprimary->Yv(),aodprimary->Zv());
841 } //pure MC kine histos
845 //__________________________________________________________________
846 Int_t AliAnaElectron::GetBtag(AliAODTrack * tr )
848 //This method uses the Displaced Vertex between electron-hadron
849 //pairs and the primary vertex to determine whether an electron is
850 //likely from a B hadron.
853 for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;
854 if (ncls1 < fITSCut) return 0;
857 //Note: 02-Sep-2009, Must use GetPosition, not XYZAtDCA
858 //Bool_t gotit = tr->XYZAtDCA(x);
859 Bool_t isNotDCA = tr->GetPosition(x);
860 if(isNotDCA) { printf("##Problem getting impact parameter!\n"); return 0; }
862 Double_t d1 = TMath::Sqrt(x[0]*x[0] + x[1]*x[1]);
863 if (TMath::Abs(d1) > fImpactCut ) return 0;
864 if (TMath::Abs(x[2]) > fImpactCut ) return 0;
865 //printf("----- impact parameter: x=%f, y=%f, z=%f -------\n",x[0],x[1], x[2]);
871 for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {
873 AliAODTrack* track2 = (AliAODTrack*) (GetAODCTS()->At(k2));
874 Int_t id1 = tr->GetID();
875 Int_t id2 = track2->GetID();
876 if(id1 == id2) continue;
879 for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
880 if (ncls2 < fITSCut) return 0;
882 if(track2->Pt() < fAssocPtCut) continue;
884 Double_t dphi = tr->Phi() - track2->Phi();
885 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
886 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
887 Double_t deta = tr->Eta() - track2->Eta();
888 Double_t dr = sqrt(deta*deta + dphi*dphi);
890 if(dr > fDrCut) continue;
892 Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);
893 if (sDca1 > fSdcaCut) nvtx1++;
894 Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);
895 if (sDca2 > fSdcaCut) nvtx2++;
896 Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);
897 if (sDca3 > fSdcaCut) nvtx3++;
899 } //loop over hadrons
902 if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);
903 if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);
904 if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);
908 fhBtagCut1->Fill(nvtx1,tr->Pt());
909 fhBtagCut2->Fill(nvtx2,tr->Pt());
910 fhBtagCut3->Fill(nvtx3,tr->Pt());
915 //__________________________________________________________________
916 Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)
918 //Compute the signed dca between two tracks
919 //and return the result
921 Double_t signDca=-999.;
922 if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);
924 //=====Now calculate DCA between both tracks=======
925 Double_t massE = 0.000511;
926 Double_t massK = 0.493677;
928 Double_t bfield = 5.; //kG
929 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
931 Double_t vertex[3] = {-999.,-999.,-999}; //vertex
932 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {
933 GetReader()->GetVertex(vertex); //If only one file, get the vertex from there
934 //FIXME: Add a check for whether file 2 is PYTHIA or HIJING
935 //If PYTHIA, then set the vertex from file 2, if not, use the
937 if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);
940 TVector3 primV(vertex[0],vertex[1],vertex[2]) ;
942 if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;
944 AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);
945 AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);
947 Double_t xplane1 = 0.; Double_t xplane2 = 0.;
948 Double_t pairdca = param1->GetDCA(param2,bfield,xplane1,xplane2);
950 Int_t id1 = 0, id2 = 0;
951 AliESDv0 bvertex(*param1,id1,*param2,id2);
953 bvertex.GetXYZ(vx,vy,vz);
957 param1->PxPyPz(emom);
958 param2->PxPyPz(hmom);
959 TVector3 emomAtB(emom[0],emom[1],emom[2]);
960 TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);
961 TVector3 secvtxpt(vx,vy,vz);
962 TVector3 decayvector(0,0,0);
963 decayvector = secvtxpt - primV; //decay vector from PrimVtx
964 Double_t decaylength = decayvector.Mag();
967 printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );
968 printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );
971 if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {
972 TVector3 sumMom = emomAtB+hmomAtB;
973 Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);
974 Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);
975 Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);
976 Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
977 Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));
978 if (mass > masscut && massPhot > 0.1) signDca = decayvector.Dot(emomAtB)/emomAtB.Mag();
979 if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);
980 if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);
990 //__________________________________________________________________
991 Bool_t AliAnaElectron::IsItPhotonic(const AliAODPWG4Particle* part)
993 //This method checks the opening angle and invariant mass of
994 //electron pairs to see if they are likely to be photonic electrons
996 Bool_t itIS = kFALSE;
998 Double_t massE = 0.000511;
999 Double_t bfield = 5.; //kG
1000 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
1002 Int_t trackId = part->GetTrackLabel(0);
1003 AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);
1004 Int_t pdg1 = part->GetPdg();
1006 AliExternalTrackParam *param1 = new AliExternalTrackParam(track);
1008 //Loop on stored AOD electrons and compute the angle differences and Minv
1009 for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {
1010 AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);
1011 Int_t pdg2 = part2->GetPdg();
1012 if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;
1013 if(part2->GetDetector() != fCalorimeter) continue;
1015 //JLK: Check opp. sign pairs only ?
1016 if(pdg1*pdg2 < 0) continue;
1018 //propagate to common vertex and check opening angle
1019 Int_t track2Id = part2->GetTrackLabel(0);
1020 AliExternalTrackParam *param2 = new AliExternalTrackParam((AliAODTrack*)GetAODCTS()->At(track2Id));
1021 Int_t id1 = 0, id2 = 0;
1022 AliESDv0 photonVtx(*param1,id1,*param2,id2);
1024 photonVtx.GetXYZ(vx,vy,vz);
1028 param1->PxPyPz(p1mom);
1029 param2->PxPyPz(p2mom);
1031 TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);
1032 TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);
1033 TVector3 sumMom = p1momAtB+p2momAtB;
1035 Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);
1036 Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);
1037 Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
1039 Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);
1040 fh1OpeningAngle->Fill(dphi);
1041 fh1MinvPhoton->Fill(mass);
1044 if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");
1058 //__________________________________________________________________
1059 void AliAnaElectron::Print(const Option_t * opt) const
1061 //Print some relevant parameters set for the analysis
1066 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1067 AliAnaPartCorrBaseClass::Print(" ");
1069 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1070 printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax);
1071 printf("residual cut = %f\n",fResidualCut);
1072 printf("---Btagging\n");
1073 printf("max IP-cut (e,h) = %f\n",fImpactCut);
1074 printf("min ITS-hits = %d\n",fITSCut);
1075 printf("max dR (e,h) = %f\n",fDrCut);
1076 printf("max pairDCA = %f\n",fPairDcaCut);
1077 printf("max decaylength = %f\n",fDecayLenCut);
1078 printf("min Associated Pt = %f\n",fAssocPtCut);
1083 //________________________________________________________________________
1084 void AliAnaElectron::ReadHistograms(TList* outputList)
1086 // Needed when Terminate is executed in distributed environment
1087 // Refill analysis histograms of this class with corresponding
1088 // histograms in output list.
1090 // Histograms of this analsys are kept in the same list as other
1091 // analysis, recover the position of
1092 // the first one and then add the next
1093 Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));
1095 //Read histograms, must be in the same order as in
1096 //GetCreateOutputObject.
1097 fh1pOverE = (TH1F *) outputList->At(index);
1098 fh1dR = (TH1F *) outputList->At(index++);
1099 fh2EledEdx = (TH2F *) outputList->At(index++);
1100 fh2MatchdEdx = (TH2F *) outputList->At(index++);
1104 //__________________________________________________________________
1105 void AliAnaElectron::Terminate(TList* outputList)
1108 //Do some plots to end
1109 //Recover histograms from output histograms list, needed for
1110 //distributed analysis.
1111 //ReadHistograms(outputList);
1113 printf(" AliAnaElectron::Terminate() *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;