]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALIsoPhoton.cxx
CommitLineData
bd092f0f 1// $Id$
30159e6f 2
bd092f0f 3#include "AliAnalysisTaskEMCALIsoPhoton.h"
30159e6f 4
bd092f0f 5#include <TCanvas.h>
6#include <TChain.h>
7#include <TFile.h>
8#include <TH1F.h>
9#include <TH2F.h>
caaf99d3 10#include <TH3F.h>
965c985f 11#include <THnSparse.h>
bd092f0f 12#include <TLorentzVector.h>
a62631a9 13#include <TList.h>
bd092f0f 14
15#include "AliAnalysisManager.h"
16#include "AliAnalysisTask.h"
17#include "AliEMCALGeometry.h"
18#include "AliEMCALRecoUtils.h"
19#include "AliESDCaloCells.h"
20#include "AliESDCaloCluster.h"
30159e6f 21#include "AliESDEvent.h"
22#include "AliESDHeader.h"
30159e6f 23#include "AliESDInputHandler.h"
bd092f0f 24#include "AliESDUtils.h"
30159e6f 25#include "AliESDtrack.h"
26#include "AliESDtrackCuts.h"
2e20efe5 27#include "AliAODEvent.h"
28#include "AliAODTrack.h"
bd092f0f 29#include "AliMCEvent.h"
30#include "AliMCEventHandler.h"
31#include "AliStack.h"
68faccf4 32#include "AliAODMCParticle.h"
2e20efe5 33#include "AliVEvent.h"
34#include "AliVTrack.h"
30159e6f 35#include "AliV0vertexer.h"
30159e6f 36#include "AliVCluster.h"
822c9944 37#include "AliOADBContainer.h"
30159e6f 38
c2aad3ae 39#include <iostream>
40using std::cout;
41using std::endl;
42
30159e6f 43ClassImp(AliAnalysisTaskEMCALIsoPhoton)
44
45//________________________________________________________________________
bd092f0f 46AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
47 AliAnalysisTaskSE(),
2e20efe5 48 fESDClusters(0),
49 fAODClusters(0),
bd092f0f 50 fSelPrimTracks(0),
3f4073ba 51 fTracks(0),
68faccf4 52 fAODMCParticles(0),
2e20efe5 53 fESDCells(0),
54 fAODCells(0),
bd092f0f 55 fPrTrCuts(0),
c45212c6 56 fCompTrCuts(0),
bd092f0f 57 fGeom(0x0),
9148fa89 58 fGeoName("EMCAL_COMPLETEV1"),
4d4ce2d2 59 fOADBContainer(0),
503a68b2 60 fVecPv(0.,0.,0.),
9148fa89 61 fPeriod("LHC11c"),
62 fTrigBit("kEMC7"),
63 fIsTrain(0),
c7bb0b43 64 fIsMc(0),
ecd47673 65 fDebug(0),
f3843637 66 fPathStrOpt("/"),
9148fa89 67 fExoticCut(0.97),
68 fIsoConeR(0.4),
965c985f 69 fNDimensions(7),
70 fECut(3.),
16a4050e 71 fTrackMult(0),
f507c09b 72 fMcIdFamily(""),
73 fNClusForDirPho(0),
74 fDirPhoPt(0),
f9e2362a 75 fHigherPtCone(0),
26a8db0f 76 fImportGeometryFromFile(0),
77 fImportGeometryFilePath(""),
2b7205ad 78 fMaxPtTrack(0),
79 fMaxEClus(0),
112eb594 80 fNCells50(0),
48aab590 81 fFilterBit(0),
98a8397b 82 fSelHybrid(kFALSE),
fe617021 83 fFillQA(kFALSE),
ed39f27f 84 fClusIdFromTracks(""),
1c764315 85 fCpvFromTrack(kFALSE),
12387380 86 fNBinsPt(200),
af4bc7bc 87 fPtBinLowEdge(0),
88 fPtBinHighEdge(100),
5bdfdbea 89 fRemMatchClus(kFALSE),
90 fMinIsoClusE(0),
877552c2 91 fNCuts(5),
d1cc86fd 92 fTrCoreRem(kFALSE),
d9cc477d 93 fClusTDiff(30e-9),
6300a3af 94 fPileUpRejSPD(kFALSE),
af4bc7bc 95 fDistToBadChan(0),
7022f12d 96 fInConeInvMass(""),
e38ddacc 97 fInConePairClEt(""),
bd092f0f 98 fESD(0),
2e20efe5 99 fAOD(0),
85b52a52 100 fVEvent(0),
c7bb0b43 101 fMCEvent(0),
102 fStack(0),
bd092f0f 103 fOutputList(0),
104 fEvtSel(0),
0b21f686 105 fNClusEt10(0),
3aca451b 106 fClusArrayNames(0),
cc57d293 107 fRecoPV(0),
f507c09b 108 fPVtxZ(0),
109 fTrMultDist(0),
1495ca60 110 fClusEtCPVSBGISO(0),
111 fClusEtCPVBGISO(0),
caaf99d3 112 fMCDirPhotonPtEtaPhi(0),
e53ab710 113 fMCIsoDirPhotonPtEtaPhi(0),
dce6be79 114 fMCDirPhotonPtEtIso(0),
c7bb0b43 115 fDecayPhotonPtMC(0),
bd092f0f 116 fCellAbsIdVsAmpl(0),
16a4050e 117 fNClusHighClusE(0),
f9e2362a 118 fHigherPtConeM02(0),
f507c09b 119 fClusEtMcPt(0),
120 fClusMcDetaDphi(0),
121 fNClusPerPho(0),
f912f9a9 122 fMcPtInConeBG(0),
123 fMcPtInConeSBG(0),
124 fMcPtInConeBGnoUE(0),
125 fMcPtInConeSBGnoUE(0),
559660d1 126 fMcPtInConeTrBGnoUE(0),
127 fMcPtInConeTrSBGnoUE(0),
1495ca60 128 fMcPtInConeMcPhoPt(0),
81f2660f 129 fAllIsoEtMcGamma(0),
130 fAllIsoNoUeEtMcGamma(0),
092ceec8 131 fMCDirPhotonPtEtaPhiNoClus(0),
e38ddacc 132 fEtCandIsoAndIsoWoPairEt(0),
4ec9dd81 133 fInConePairedClusEtVsCandEt(0),
2b7205ad 134 fHnOutput(0),
135 fQAList(0),
136 fNTracks(0),
137 fEmcNCells(0),
138 fEmcNClus(0),
139 fEmcNClusCut(0),
140 fNTracksECut(0),
141 fEmcNCellsCut(0),
ed39f27f 142 fEmcClusETM1(0),
143 fEmcClusETM2(0),
34393c41 144 fEmcClusNotExo(0),
877552c2 145 fEmcClusEClusCuts(0),
e4ba80ad 146 fEmcClusEPhi(0),
147 fEmcClusEPhiCut(0),
148 fEmcClusEEta(0),
149 fEmcClusEEtaCut(0),
150 fTrackPtPhi(0),
151 fTrackPtPhiCut(0),
152 fTrackPtEta(0),
3ae97198 153 fTrackPtEtaCut(0),
74d5d8ab 154 fMaxCellEPhi(0),
ed4488b6 155 fDetaDphiFromTM(0),
156 fEoverPvsE(0)
bd092f0f 157{
158 // Default constructor.
26a8db0f 159 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
bd092f0f 160}
161
162//________________________________________________________________________
163AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
164 AliAnalysisTaskSE(name),
2e20efe5 165 fESDClusters(0),
166 fAODClusters(0),
30159e6f 167 fSelPrimTracks(0),
3f4073ba 168 fTracks(0),
68faccf4 169 fAODMCParticles(0),
2e20efe5 170 fESDCells(0),
171 fAODCells(0),
30159e6f 172 fPrTrCuts(0),
c45212c6 173 fCompTrCuts(0),
30159e6f 174 fGeom(0x0),
175 fGeoName("EMCAL_COMPLETEV1"),
4d4ce2d2 176 fOADBContainer(0),
99f97ac0 177 fVecPv(0.,0.,0.),
30159e6f 178 fPeriod("LHC11c"),
751194e8 179 fTrigBit("kEMC7"),
30159e6f 180 fIsTrain(0),
c7bb0b43 181 fIsMc(0),
ecd47673 182 fDebug(0),
f3843637 183 fPathStrOpt("/"),
30159e6f 184 fExoticCut(0.97),
185 fIsoConeR(0.4),
965c985f 186 fNDimensions(7),
187 fECut(3.),
16a4050e 188 fTrackMult(0),
f507c09b 189 fMcIdFamily(""),
190 fNClusForDirPho(0),
191 fDirPhoPt(0),
f9e2362a 192 fHigherPtCone(0),
26a8db0f 193 fImportGeometryFromFile(0),
194 fImportGeometryFilePath(""),
2b7205ad 195 fMaxPtTrack(0),
196 fMaxEClus(0),
112eb594 197 fNCells50(0),
48aab590 198 fFilterBit(0),
98a8397b 199 fSelHybrid(kFALSE),
fe617021 200 fFillQA(kFALSE),
ed39f27f 201 fClusIdFromTracks(""),
1c764315 202 fCpvFromTrack(kFALSE),
12387380 203 fNBinsPt(200),
af4bc7bc 204 fPtBinLowEdge(0.),
205 fPtBinHighEdge(100),
5bdfdbea 206 fRemMatchClus(kFALSE),
207 fMinIsoClusE(0),
877552c2 208 fNCuts(5),
d1cc86fd 209 fTrCoreRem(kFALSE),
d9cc477d 210 fClusTDiff(30e-9),
6300a3af 211 fPileUpRejSPD(kFALSE),
af4bc7bc 212 fDistToBadChan(0),
7022f12d 213 fInConeInvMass(""),
e38ddacc 214 fInConePairClEt(""),
30159e6f 215 fESD(0),
2e20efe5 216 fAOD(0),
85b52a52 217 fVEvent(0),
c7bb0b43 218 fMCEvent(0),
219 fStack(0),
30159e6f 220 fOutputList(0),
30159e6f 221 fEvtSel(0),
0b21f686 222 fNClusEt10(0),
3aca451b 223 fClusArrayNames(0),
cc57d293 224 fRecoPV(0),
bd092f0f 225 fPVtxZ(0),
f507c09b 226 fTrMultDist(0),
1495ca60 227 fClusEtCPVSBGISO(0),
228 fClusEtCPVBGISO(0),
caaf99d3 229 fMCDirPhotonPtEtaPhi(0),
e53ab710 230 fMCIsoDirPhotonPtEtaPhi(0),
dce6be79 231 fMCDirPhotonPtEtIso(0),
c7bb0b43 232 fDecayPhotonPtMC(0),
233 fCellAbsIdVsAmpl(0),
bd092f0f 234 fNClusHighClusE(0),
f9e2362a 235 fHigherPtConeM02(0),
f507c09b 236 fClusEtMcPt(0),
237 fClusMcDetaDphi(0),
238 fNClusPerPho(0),
f912f9a9 239 fMcPtInConeBG(0),
240 fMcPtInConeSBG(0),
241 fMcPtInConeBGnoUE(0),
242 fMcPtInConeSBGnoUE(0),
559660d1 243 fMcPtInConeTrBGnoUE(0),
244 fMcPtInConeTrSBGnoUE(0),
1495ca60 245 fMcPtInConeMcPhoPt(0),
81f2660f 246 fAllIsoEtMcGamma(0),
247 fAllIsoNoUeEtMcGamma(0),
092ceec8 248 fMCDirPhotonPtEtaPhiNoClus(0),
e38ddacc 249 fEtCandIsoAndIsoWoPairEt(0),
4ec9dd81 250 fInConePairedClusEtVsCandEt(0),
2b7205ad 251 fHnOutput(0),
252 fQAList(0),
253 fNTracks(0),
254 fEmcNCells(0),
255 fEmcNClus(0),
256 fEmcNClusCut(0),
257 fNTracksECut(0),
258 fEmcNCellsCut(0),
ed39f27f 259 fEmcClusETM1(0),
260 fEmcClusETM2(0),
34393c41 261 fEmcClusNotExo(0),
877552c2 262 fEmcClusEClusCuts(0),
e4ba80ad 263 fEmcClusEPhi(0),
264 fEmcClusEPhiCut(0),
265 fEmcClusEEta(0),
266 fEmcClusEEtaCut(0),
267 fTrackPtPhi(0),
268 fTrackPtPhiCut(0),
269 fTrackPtEta(0),
3ae97198 270 fTrackPtEtaCut(0),
74d5d8ab 271 fMaxCellEPhi(0),
ed4488b6 272 fDetaDphiFromTM(0),
273 fEoverPvsE(0)
30159e6f 274{
275 // Constructor
276
277 // Define input and output slots here
278 // Input slot #0 works with a TChain
279 DefineInput(0, TChain::Class());
280 // Output slot #0 id reserved by the base class for AOD
281 // Output slot #1 writes into a TH1 container
282 DefineOutput(1, TList::Class());
2b7205ad 283 DefineOutput(2, TList::Class());
30159e6f 284}
bd092f0f 285
30159e6f 286//________________________________________________________________________
287void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
288{
bd092f0f 289 // Create histograms, called once.
30159e6f 290
2e20efe5 291 fESDClusters = new TObjArray();
652a3f8f 292 fAODClusters = new TObjArray();
30159e6f 293 fSelPrimTracks = new TObjArray();
294
295
296 fOutputList = new TList();
a62631a9 297 fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
30159e6f 298
f507c09b 299 fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
4d4ce2d2 300 fOADBContainer = new AliOADBContainer("AliEMCALgeo");
301 fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
302
30159e6f 303 fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
304 fOutputList->Add(fEvtSel);
0b21f686 305
306 fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
307 fOutputList->Add(fNClusEt10);
3aca451b 308
309 fClusArrayNames = new TH1F("hClusArrayNames","cluster array names (0=CaloClusters,1=EmcCaloClusters,2=Others);option;#events",3,0,3);
310 fOutputList->Add(fClusArrayNames);
30159e6f 311
cc57d293 312 fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
313 fOutputList->Add(fRecoPV);
314
30159e6f 315 fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
316 fOutputList->Add(fPVtxZ);
c7bb0b43 317
f507c09b 318 fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
319 fOutputList->Add(fTrMultDist);
320
1495ca60 321 fClusEtCPVSBGISO = new TH2F("hClusEtCPVSBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.1<#lambda_{0}^{2}<0.3;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
322 fClusEtCPVSBGISO->Sumw2();
323 fOutputList->Add(fClusEtCPVSBGISO);
324
325 fClusEtCPVBGISO = new TH2F("hClusEtCPVBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.5<#lambda_{0}^{2}<2.0;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
326 fClusEtCPVBGISO->Sumw2();
327 fOutputList->Add(fClusEtCPVBGISO);
328
f4e24fd3 329 fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
caaf99d3 330 fMCDirPhotonPtEtaPhi->Sumw2();
331 fOutputList->Add(fMCDirPhotonPtEtaPhi);
c7bb0b43 332
f4e24fd3 333 fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
e53ab710 334 fMCIsoDirPhotonPtEtaPhi->Sumw2();
335 fOutputList->Add(fMCIsoDirPhotonPtEtaPhi);
336
f4e24fd3 337 fMCDirPhotonPtEtIso = new TH2F("hMCDirPhotonPtEtIso",Form("photon (gq->#gammaq @MC) p_{T}, E_{T}^{ISO} (R=%1.1f);GeV/c;E_{T}^{ISO} GeV/c",fIsoConeR),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,20,-0.25,9.75);
dce6be79 338 fMCDirPhotonPtEtIso->Sumw2();
339 fOutputList->Add(fMCDirPhotonPtEtIso);
340
341
12387380 342 fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
c7bb0b43 343 fDecayPhotonPtMC->Sumw2();
344 fOutputList->Add(fDecayPhotonPtMC);
345
30159e6f 346 fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
a62631a9 347 fOutputList->Add(fCellAbsIdVsAmpl);
30159e6f 348
349 fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
350 fOutputList->Add(fNClusHighClusE);
351
12387380 352 fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,400,0,4);
f9e2362a 353 fOutputList->Add(fHigherPtConeM02);
354
f507c09b 355 fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
356 fOutputList->Add(fClusEtMcPt);
357
358 fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
359 fOutputList->Add(fClusMcDetaDphi);
360
361 fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
362 fOutputList->Add(fNClusPerPho);
363
d4f449df 364 fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
f912f9a9 365 fOutputList->Add(fMcPtInConeBG);
366
d4f449df 367 fMcPtInConeSBG = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
f912f9a9 368 fOutputList->Add(fMcPtInConeSBG);
369
d4f449df 370 fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
f912f9a9 371 fOutputList->Add(fMcPtInConeBGnoUE);
372
d4f449df 373 fMcPtInConeSBGnoUE = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
f912f9a9 374 fOutputList->Add(fMcPtInConeSBGnoUE);
375
559660d1 376 fMcPtInConeTrBGnoUE = new TH2F("hMcPtInConeTrBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (BG template);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
377 fOutputList->Add(fMcPtInConeTrBGnoUE);
378
379 fMcPtInConeTrSBGnoUE = new TH2F("hMcPtInConeTrSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (SBG range);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
380 fOutputList->Add(fMcPtInConeTrSBGnoUE);
381
1495ca60 382 fMcPtInConeMcPhoPt = new TH2F("hMcPtInConeMcPhoPt","#sum_{in-cone}p_{T}^{mc-primaries} vs. prompt photon p_{T};p_{T}^{mc-#gamma} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,-0.25,99.75);
383 fOutputList->Add(fMcPtInConeMcPhoPt);
384
12387380 385 fAllIsoEtMcGamma = new TH2F("hAllIsoEtMcGammaE","ISO^{TRK+EMC} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
81f2660f 386 fOutputList->Add(fAllIsoEtMcGamma);
387
12387380 388 fAllIsoNoUeEtMcGamma = new TH2F("hAllIsoNoUeEtMcGammaE","ISO^{TRK+EMC}_{noue} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC}_{noue} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
81f2660f 389 fOutputList->Add(fAllIsoNoUeEtMcGamma);
390
f912f9a9 391
12387380 392 fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
092ceec8 393 fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
f507c09b 394
e38ddacc 395 fEtCandIsoAndIsoWoPairEt = new TH3F("hEtCandIsoAndIsoWoPairEt","E_{T}^{cand} vs. E_{T}^{ISO} (EMC+Trk) (0.1<M02<0.3, 0.110<m_{#gamma#gamma}<0.165 only);E_{T}^{cand}; E_{T}^{ISO}; E_{T}^{ISO} (w/o #pi^{0} pair E_{T})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,200,1000,0,200);
396 fOutputList->Add(fEtCandIsoAndIsoWoPairEt);
7022f12d 397
4ec9dd81 398 fInConePairedClusEtVsCandEt = new TH2F("hInConePairedClusEtVsCandEt","E_{T}^{partner} vs. E_{T}^{cand} (R<0.4, 0.110<m_{#gamma#gamma}<0.165);E_{T}^{cand};E_{T}^{partner}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,0,40);
399 fOutputList->Add(fInConePairedClusEtVsCandEt);
400
c141c3fd 401 Int_t nEt=fNBinsPt*5, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nCeIsoNoUE=1000, nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=fNBinsPt, nInConeMass=100;
402 Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt,nInConeMass};
c1f18270 403 fNDimensions = sizeof(bins)/sizeof(Int_t);
404 const Int_t ndims = fNDimensions;
c141c3fd 405 Double_t xmin[] = { fPtBinLowEdge, 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,fPtBinLowEdge,0.0};
406 Double_t xmax[] = { fPtBinHighEdge, 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,fPtBinHighEdge, 1.0};
80c98845 407 if(fPeriod.Contains("11h")){
408 xmax[12]=3999.5;
409 }
f507c09b 410 fHnOutput = new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
c141c3fd 411 fHnOutput->Sumw2();
965c985f 412 fOutputList->Add(fHnOutput);
413
2b7205ad 414 //QA outputs
415 fQAList = new TList();
416 fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
965c985f 417
2b7205ad 418 fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5);
419 fNTracks->Sumw2();
420 fQAList->Add(fNTracks);
421
422 fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5);
423 fEmcNCells->Sumw2();
424 fQAList->Add(fEmcNCells);
425 fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5);
426 fEmcNClus->Sumw2();
427 fQAList->Add(fEmcNClus);
2074fcf3 428 fEmcNClusCut = new TH1F("fEmcNClusCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
2b7205ad 429 fEmcNClusCut->Sumw2();
430 fQAList->Add(fEmcNClusCut);
2074fcf3 431 fNTracksECut = new TH1F("fNTracksECut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
2b7205ad 432 fNTracksECut->Sumw2();
433 fQAList->Add(fNTracksECut);
2074fcf3 434 fEmcNCellsCut = new TH1F("fEmcNCellsCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
2b7205ad 435 fEmcNCellsCut->Sumw2();
436 fQAList->Add(fEmcNCellsCut);
12387380 437 fEmcClusETM1 = new TH1F("fEmcClusETM1","(method clus->GetTrackDx,z);GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
ed39f27f 438 fEmcClusETM1->Sumw2();
439 fQAList->Add(fEmcClusETM1);
12387380 440 fEmcClusETM2 = new TH1F("fEmcClusETM2","(method track->GetEMCALcluster());GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
ed39f27f 441 fEmcClusETM2->Sumw2();
442 fQAList->Add(fEmcClusETM2);
34393c41 443 fEmcClusNotExo = new TH1F("fEmcClusNotExo","exotics removed;GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
444 fEmcClusNotExo->Sumw2();
445 fQAList->Add(fEmcClusNotExo);
12387380 446 fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
e4ba80ad 447 fEmcClusEPhi->Sumw2();
448 fQAList->Add(fEmcClusEPhi);
2074fcf3 449 fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",Form("(at least one E_{clus}>%1.1f);GeV;#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
e4ba80ad 450 fEmcClusEPhiCut->Sumw2();
451 fQAList->Add(fEmcClusEPhiCut);
12387380 452 fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,19,-0.9,0.9);
e4ba80ad 453 fEmcClusEEta->Sumw2();
454 fQAList->Add(fEmcClusEEta);
2074fcf3 455 fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",Form("(at least one E_{clus}>%1.1f);GeV;#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
e4ba80ad 456 fEmcClusEEtaCut->Sumw2();
457 fQAList->Add(fEmcClusEEtaCut);
12387380 458 fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
e4ba80ad 459 fTrackPtPhi->Sumw2();
460 fQAList->Add(fTrackPtPhi);
2074fcf3 461 fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
e4ba80ad 462 fTrackPtPhiCut->Sumw2();
463 fQAList->Add(fTrackPtPhiCut);
12387380 464 fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
e4ba80ad 465 fTrackPtEta->Sumw2();
466 fQAList->Add(fTrackPtEta);
2074fcf3 467 fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
e4ba80ad 468 fTrackPtEtaCut->Sumw2();
469 fQAList->Add(fTrackPtEtaCut);
877552c2 470 fEmcClusEClusCuts = new TH2F("fEmcClusEClusCuts",";GeV;cut",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,fNCuts,-0.5,fNCuts-0.5);
471 fEmcClusEClusCuts->Sumw2();
472 fQAList->Add(fEmcClusEClusCuts);
3ae97198 473
12387380 474 fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
3ae97198 475 fMaxCellEPhi->Sumw2();
476 fQAList->Add(fMaxCellEPhi);
477
74d5d8ab 478 fDetaDphiFromTM = new TH2F("fDetaDphiFromTM","d#phi vs. d#eta of clusters from track->GetEMCALcluster();d#eta;d#phi",100,-0.05,0.05,200,-0.1,0.1);
479 fDetaDphiFromTM->Sumw2();
480 fQAList->Add(fDetaDphiFromTM);
481
ed4488b6 482 fEoverPvsE = new TH2F("fEoverPvsE","E^{clus}/p^{track} vs E^{clus} (80<TPCsignal<100);E^{clus} [GeV];E^{clus}/p^{track} [c^{-1}]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,100,0,2);
483 fEoverPvsE->Sumw2();
484 fQAList->Add(fEoverPvsE);
485
30159e6f 486 PostData(1, fOutputList);
2b7205ad 487 PostData(2, fQAList);
30159e6f 488}
489
490//________________________________________________________________________
491void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
492{
bd092f0f 493 // Main loop, called for each event.
2e20efe5 494 fESDClusters = 0;
495 fESDCells = 0;
496 fAODClusters = 0;
497 fAODCells = 0;
bd092f0f 498 // event trigger selection
30159e6f 499 Bool_t isSelected = 0;
751194e8 500 if(fPeriod.Contains("11a")){
501 if(fTrigBit.Contains("kEMC"))
502 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
503 else
504 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
505 }
506 else{
507 if(fTrigBit.Contains("kEMC"))
508 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
509 else
f507c09b 510 if(fTrigBit.Contains("kMB"))
511 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
512 else
513 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
751194e8 514 }
80c98845 515 if(fPeriod.Contains("11h")){
516 if(fTrigBit.Contains("kAny"))
517 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
518 if(fTrigBit.Contains("kAnyINT"))
519 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
520 }
c7bb0b43 521 if(fIsMc)
522 isSelected = kTRUE;
ecd47673 523 if(fDebug)
524 printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
30159e6f 525 if(!isSelected )
526 return;
f3843637 527 if(fIsMc){
528 TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
529 TFile *file = (TFile*)tree->GetCurrentFile();
530 TString filename = file->GetName();
531 if(!filename.Contains(fPathStrOpt.Data()))
532 return;
533 }
85b52a52 534 fVEvent = (AliVEvent*)InputEvent();
535 if (!fVEvent) {
2e20efe5 536 printf("ERROR: event not available\n");
30159e6f 537 return;
538 }
822c9944 539 Int_t runnumber = InputEvent()->GetRunNumber() ;
7a7744db 540 if(fDebug)
541 printf("run number = %d\n",runnumber);
542
85b52a52 543 fESD = dynamic_cast<AliESDEvent*>(fVEvent);
7a7744db 544 if(!fESD){
85b52a52 545 fAOD = dynamic_cast<AliAODEvent*>(fVEvent);
7a7744db 546 if(!fAOD){
547 printf("ERROR: Invalid type of event!!!\n");
548 return;
549 }
550 else if(fDebug)
551 printf("AOD EVENT!\n");
552 }
30159e6f 553
554 fEvtSel->Fill(0);
ecd47673 555 if(fDebug)
26a8db0f 556 printf("event is ok,\n run number=%d\n",runnumber);
557
30159e6f 558
85b52a52 559 AliVVertex *pv = (AliVVertex*)fVEvent->GetPrimaryVertex();
2e20efe5 560 Bool_t pvStatus = kTRUE;
561 if(fESD){
562 AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
563 pvStatus = esdv->GetStatus();
564 }
425adc2f 565 /*if(fAOD){
566 AliAODVertex *aodv = (AliAODVertex*)fAOD->GetPrimaryVertex();
567 pvStatus = aodv->GetStatus();
568 }*/
d1582697 569 if(!pv)
30159e6f 570 return;
2e20efe5 571 if(!pvStatus)
d1582697 572 fRecoPV->Fill(0);
cc57d293 573 else
574 fRecoPV->Fill(1);
30159e6f 575 fPVtxZ->Fill(pv->GetZ());
503a68b2 576 fVecPv.SetXYZ(pv->GetX(),pv->GetY(),pv->GetZ());
425adc2f 577 if(TMath::Abs(pv->GetZ())>10)
30159e6f 578 return;
ecd47673 579 if(fDebug)
580 printf("passed vertex cut\n");
30159e6f 581
582 fEvtSel->Fill(1);
6300a3af 583 if(fVEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.) && fPileUpRejSPD){
584 if(fDebug)
585 printf("Event is SPD pile-up!***\n");
586 return;
587 }
2e20efe5 588 if(fESD)
589 fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
590 if(fAOD)
591 fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
bd092f0f 592
06a28959 593 if(!fTracks){
594 AliError("Track array in event is NULL!");
2e20efe5 595 if(fDebug)
596 printf("returning due to not finding tracks!\n");
06a28959 597 return;
598 }
30159e6f 599 // Track loop to fill a pT spectrum
3f4073ba 600 const Int_t Ntracks = fTracks->GetEntriesFast();
601 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
602 // for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
603 //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
2e20efe5 604 AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
30159e6f 605 if (!track)
606 continue;
2e20efe5 607 AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
608 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
ed39f27f 609 if(esdTrack){
610 if(esdTrack->GetEMCALcluster()>0)
611 fClusIdFromTracks.Append(Form("%d ",esdTrack->GetEMCALcluster()));
612 if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
613 fSelPrimTracks->Add(track);
c45212c6 614 } else if(fCompTrCuts && fCompTrCuts->IsSelected(track)) {
615 fSelPrimTracks->Add(track);
ed39f27f 616 }
30159e6f 617 }
2b7205ad 618 else if(aodTrack){
98a8397b 619 if (fSelHybrid && !aodTrack->IsHybridGlobalConstrainedGlobal())
dcd3dca1 620 continue ;
98a8397b 621 if(!fSelHybrid && !aodTrack->TestFilterBit(fFilterBit))
622 continue;
2e20efe5 623 fSelPrimTracks->Add(track);
2b7205ad 624 /*if(fTrackMaxPt<track->Pt())
625 fTrackMaxPt = track->Pt();*/
626 }
bd092f0f 627 }
628
4d4ce2d2 629 TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
822c9944 630 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
631 if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
632 break;
633 /*if(fESD)
634 fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
635 else*/
85b52a52 636 // if(fVEvent->GetEMCALMatrix(mod))
822c9944 637 fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
638 fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
30159e6f 639 }
4d4ce2d2 640
2e20efe5 641 if(fESD){
642 AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
643 fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
644 if(fDebug)
645 printf("ESD Track mult= %d\n",fTrackMult);
646 }
647 else if(fAOD){
648 fTrackMult = Ntracks;
649 if(fDebug)
650 printf("AOD Track mult= %d\n",fTrackMult);
651 }
f507c09b 652 fTrMultDist->Fill(fTrackMult);
3aca451b 653 TList *l = 0;
654 TString clusArrayName = "";
2e20efe5 655 if(fESD){
3aca451b 656 l = fESD->GetList();
731bd9d1 657 if(fDebug)
658 l->Print();
3aca451b 659 for(int nk=0;nk<l->GetEntries();nk++){
660 TObject *obj = (TObject*)l->At(nk);
661 TString oname = obj->GetName();
662 if(oname.Contains("CaloClus"))
663 clusArrayName = oname;
731bd9d1 664 else
665 continue;
666 if(clusArrayName=="CaloClusters")
667 fClusArrayNames->Fill(0);
668 else{
669 if(clusArrayName=="EmcCaloClusters")
670 fClusArrayNames->Fill(1);
671 else
672 fClusArrayNames->Fill(2);
673 }
950b5606 674 }
3aca451b 675 fESDClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
2e20efe5 676 fESDCells = fESD->GetEMCALCells();
677 if(fDebug)
678 printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
679 }
680 else if(fAOD){
3aca451b 681 l = fAOD->GetList();
731bd9d1 682 if(fDebug)
683 l->Print();
3aca451b 684 //fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
685 for(int nk=0;nk<l->GetEntries();nk++){
686 TObject *obj = (TObject*)l->At(nk);
687 TString oname = obj->GetName();
688 if(oname.Contains("aloClus"))
689 clusArrayName = oname;
731bd9d1 690 else
691 continue;
692 if(clusArrayName=="caloClusters")
693 fClusArrayNames->Fill(0);
694 else{
695 if(clusArrayName=="EmcCaloClusters")
696 fClusArrayNames->Fill(1);
697 else
698 fClusArrayNames->Fill(2);
699 }
3aca451b 700 }
701 fAODClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
2e20efe5 702 fAODCells = fAOD->GetEMCALCells();
703 if(fDebug)
704 printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
705 }
3aca451b 706 if(fDebug){
731bd9d1 707 printf("clus array is named %s +++++++++\n",clusArrayName.Data());
3aca451b 708 }
731bd9d1 709
710
c7bb0b43 711 fMCEvent = MCEvent();
ecd47673 712 if(fMCEvent){
713 if(fDebug)
ec9ab86b 714 std::cout<<"MCevent exists\n";
c7bb0b43 715 fStack = (AliStack*)fMCEvent->Stack();
68faccf4 716 if(!fStack)
717 fAODMCParticles = (TClonesArray*)fVEvent->FindListObject("mcparticles");
ecd47673 718 }
719 else{
720 if(fDebug)
ec9ab86b 721 std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
ecd47673 722 }
112eb594 723 LoopOnCells();
f507c09b 724 FollowGamma();
725 if(fDebug)
726 printf("passed calling of FollowGamma\n");
727 FillClusHists();
728 if(fDebug)
729 printf("passed calling of FillClusHists\n");
63e40cb8 730 FillMcHists();
ecd47673 731 if(fDebug)
732 printf("passed calling of FillMcHists\n");
fe617021 733 if(fFillQA)
ab031886 734 FillQA();
2b7205ad 735 if(fDebug)
736 printf("passed calling of FillQA\n");
652a3f8f 737 if(fESD)
738 fESDClusters->Clear();
f507c09b 739 fSelPrimTracks->Clear();
740 fNClusForDirPho = 0;
112eb594 741 fNCells50 = 0;
ed39f27f 742 fClusIdFromTracks = "";
503a68b2 743 fVecPv.Clear();
c7bb0b43 744
30159e6f 745 PostData(1, fOutputList);
2b7205ad 746 PostData(2, fQAList);
30159e6f 747}
bd092f0f 748
30159e6f 749//________________________________________________________________________
750void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
751{
2e20efe5 752 if(fDebug)
753 printf("Inside FillClusHists()....\n");
bd092f0f 754 // Fill cluster histograms.
2e20efe5 755 TObjArray *clusters = fESDClusters;
bd092f0f 756
2e20efe5 757 if (!clusters){
758 clusters = fAODClusters;
759 if(fDebug)
760 printf("ESD clusters empty...");
761 }
762 if (!clusters){
763 if(fDebug)
764 printf(" and AOD clusters as well!!!\n");
30159e6f 765 return;
2e20efe5 766 }
767 if(fDebug)
768 printf("\n");
769
770 const Int_t nclus = clusters->GetEntries();
30159e6f 771 if(nclus==0)
772 return;
ecd47673 773 if(fDebug)
774 printf("Inside FillClusHists and there are %d clusters\n",nclus);
e681d9ce 775 Double_t maxE = 0;
0b21f686 776 Int_t nclus10 = 0;
f507c09b 777 Double_t ptmc=-1;
30159e6f 778 for(Int_t ic=0;ic<nclus;ic++){
779 maxE=0;
2e20efe5 780 AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
950b5606 781 if(!c){
782 if(fDebug)
783 printf("cluster pointer does not exist! xxxx\n");
30159e6f 784 continue;
950b5606 785 }
786 if(!c->IsEMCAL()){
787 if(fDebug)
788 printf("cluster is not EMCAL! xxxx\n");
30159e6f 789 continue;
950b5606 790 }
791 if(c->E()<fECut){
792 if(fDebug)
793 printf("cluster has E<%1.1f! xxxx\n", fECut);
965c985f 794 continue;
950b5606 795 }
796 if(fCpvFromTrack && fClusIdFromTracks.Contains(Form("%d",ic))){
797 if(fDebug)
798 printf("cluster does not pass CPV criterion! xxxx\n");
1c764315 799 continue;
950b5606 800 }
801 if(IsExotic(c)){
802 if(fDebug)
803 printf("cluster is exotic! xxxx\n");
34393c41 804 continue;
950b5606 805 }
af4bc7bc 806 if(c->GetDistanceToBadChannel()<fDistToBadChan){
807 if(fDebug)
808 printf("cluster distance to bad channel is %1.1f (<%1.1f) xxxx\n",c->GetDistanceToBadChannel(),fDistToBadChan);
809 continue;
810 }
950b5606 811 Short_t id;
30159e6f 812 Double_t Emax = GetMaxCellEnergy( c, id);
9045b7cc 813 if(fDebug)
814 printf("cluster max cell E=%1.1f",Emax);
30159e6f 815 Float_t clsPos[3] = {0,0,0};
816 c->GetPosition(clsPos);
817 TVector3 clsVec(clsPos);
9ccace04 818 clsVec -= fVecPv;
3dc2825e 819 Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
f507c09b 820 if(fDebug)
821 printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
0b21f686 822 if(Et>10)
823 nclus10++;
adce44d4 824 Float_t ceiso=0, cephiband=0, cecore=0;
825 Float_t triso=0, trphiband=0, trcore=0;
826 Float_t alliso=0, allphiband=0;//, allcore;
30159e6f 827 Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
828 Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
7022f12d 829 Bool_t isCPV = kFALSE;
830 if(TMath::Abs(c->GetTrackDx())>0.03 || TMath::Abs(c->GetTrackDz())>0.02)
831 isCPV = kTRUE;
832 GetCeIso(clsVec, id, ceiso, cephiband, cecore, Et);
30159e6f 833 GetTrIso(clsVec, triso, trphiband, trcore);
c141c3fd 834 Int_t nInConePairs = 0;
835 Double_t onePairMass = 0;
77ef8e9f 836 //---
837 //if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV){
838 TObjArray *inConeInvMassArr = (TObjArray*)fInConeInvMass.Tokenize(";");
839 TObjArray *inConePairClEt = (TObjArray*)fInConePairClEt.Tokenize(";");
840 nInConePairs = inConeInvMassArr->GetEntriesFast();
841 Int_t nInConePi0 = inConePairClEt->GetEntriesFast();
842 if(nInConePairs != nInConePi0)
843 printf("Inconsistent number of in cone pairs!!!\n");
844 for(int ipair=0;ipair<nInConePairs;ipair++){
845 TObjString *obs = (TObjString*)inConeInvMassArr->At(ipair);
846 TObjString *obet = (TObjString*)inConePairClEt->At(ipair);
847 TString smass = obs->GetString();
848 TString spairEt = obet->GetString();
849 Double_t pairmass = smass.Atof();
850 Double_t pairEt = spairEt.Atof();//this must be zero when inv mass outside pi0 range
851 if(0==ipair && nInConePairs==1)
852 onePairMass = pairmass;
853 if(fDebug)
854 printf("=================+++++++++++++++Inv mass inside the cone for photon range: %1.1f,%1.1f,%1.1f+-++++-+-+-+-++-+-+-\n",Et,pairmass,ceiso+triso);
855 fEtCandIsoAndIsoWoPairEt->Fill(Et,ceiso+triso,ceiso+triso-pairEt);
7022f12d 856 }
77ef8e9f 857 //}
858 //---
f9e2362a 859 Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
860 if(Et>10 && Et<15 && dr>0.025){
861 fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
862 if(fDebug)
863 printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
864 }
30159e6f 865 alliso = ceiso + triso;
866 allphiband = cephiband + trphiband;
6f00a63d 867 //allcore = cecore + trcore;
30159e6f 868 Float_t ceisoue = cephiband/phibandArea*netConeArea;
869 Float_t trisoue = trphiband/phibandArea*netConeArea;
870 Float_t allisoue = allphiband/phibandArea*netConeArea;
f912f9a9 871 Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
872 if(fDebug && Et>10)
873 printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
874 if(c->GetM02()>0.5 && c->GetM02()<2.0){
559660d1 875 fMcPtInConeBG->Fill(alliso-allisoue, mcptsum);
876 fMcPtInConeBGnoUE->Fill(alliso, mcptsum);
877 fMcPtInConeTrBGnoUE->Fill(triso, mcptsum);
f912f9a9 878 }
b16ded35 879 if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
559660d1 880 fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum);
881 fMcPtInConeSBGnoUE->Fill(alliso, mcptsum);
882 fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum);
81f2660f 883 if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
559660d1 884 fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue);
885 fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore);
81f2660f 886 }
f912f9a9 887 }
1495ca60 888 if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV)
d834227a 889 fClusEtCPVSBGISO->Fill(Et,alliso - trcore);
1495ca60 890 if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV)
26d2a019 891 fClusEtCPVBGISO->Fill(Et,alliso - trcore);
965c985f 892 const Int_t ndims = fNDimensions;
893 Double_t outputValues[ndims];
e53ab710 894 if(mcptsum<2)
895 ptmc = GetClusSource(c);
896 else
897 ptmc = -0.1;
965c985f 898 outputValues[0] = Et;
899 outputValues[1] = c->GetM02();
1c86c72c 900 outputValues[2] = ceiso/*cecore*/-ceisoue;
965c985f 901 outputValues[3] = triso-trisoue;
1c86c72c 902 outputValues[4] = alliso/*cecore*/-allisoue - trcore;
903 outputValues[5] = ceiso;
904 outputValues[6] = alliso - trcore;
8e574ec0 905 if(fDebug)
906 printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz());
907 if(TMath::Abs(c->GetTrackDx())<0.1)
908 outputValues[7] = c->GetTrackDx();
909 else
910 outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx());
911 if(TMath::Abs(c->GetTrackDz())<0.05)
912 outputValues[8] = c->GetTrackDz();
913 else
914 outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz());
717bb45b 915 outputValues[9] = clsVec.Eta();
916 outputValues[10] = clsVec.Phi();
2e20efe5 917 if(fESDCells)
918 outputValues[11] = fESDCells->GetCellTime(id);
919 else if(fAODCells)
920 outputValues[11] = fAODCells->GetCellTime(id);
16a4050e 921 outputValues[12] = fTrackMult;
f507c09b 922 outputValues[13] = ptmc;
c141c3fd 923 if(nInConePairs == 1)
924 outputValues[14] = onePairMass;
925 else
926 outputValues[14] = -1;
965c985f 927 fHnOutput->Fill(outputValues);
30159e6f 928 if(c->E()>maxE)
929 maxE = c->E();
930 }
931 fNClusHighClusE->Fill(maxE,nclus);
ee3d9c10 932 fMaxEClus = maxE;
0b21f686 933 fNClusEt10->Fill(nclus10);
f507c09b 934 fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
30159e6f 935}
bd092f0f 936
30159e6f 937//________________________________________________________________________
7022f12d 938void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core, Double_t EtCl)
30159e6f 939{
950b5606 940 if(fDebug)
941 printf("....indside GetCeIso funtcion\n");
bd092f0f 942 // Get cell isolation.
2e20efe5 943 AliVCaloCells *cells = fESDCells;
1c86c72c 944 if (!cells){
2e20efe5 945 cells = fAODCells;
1c86c72c 946 if(fDebug)
3aca451b 947 printf("ESD cells empty...\n");
1c86c72c 948 }
949 if (!cells){
950 if(fDebug)
951 printf(" and AOD cells are empty as well!!!\n");
30159e6f 952 return;
1c86c72c 953 }
2e20efe5 954
1c86c72c 955 TObjArray *clusters = fESDClusters;
956 if (!clusters)
957 clusters = fAODClusters;
958 if (!clusters)
959 return;
960
7022f12d 961 fInConeInvMass = "";
e38ddacc 962 fInConePairClEt="";
1c86c72c 963 const Int_t nclus = clusters->GetEntries();
964 //const Int_t ncells = cells->GetNumberOfCells();
30159e6f 965 Float_t totiso=0;
966 Float_t totband=0;
967 Float_t totcore=0;
968 Float_t etacl = vec.Eta();
969 Float_t phicl = vec.Phi();
30159e6f 970 if(phicl<0)
971 phicl+=TMath::TwoPi();
1c86c72c 972 /*Int_t absid = -1;
30159e6f 973 Float_t eta=-1, phi=-1;
974 for(int icell=0;icell<ncells;icell++){
2e20efe5 975 absid = TMath::Abs(cells->GetCellNumber(icell));
976 Float_t celltime = cells->GetCellTime(absid);
f224025f 977 //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
978 if(TMath::Abs(celltime-maxtcl)>2e-8 )
c3f33463 979 continue;
30159e6f 980 if(!fGeom)
981 return;
982 fGeom->EtaPhiFromIndex(absid,eta,phi);
983 Float_t dphi = TMath::Abs(phi-phicl);
984 Float_t deta = TMath::Abs(eta-etacl);
1c86c72c 985 Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
986 for(int ic=0;ic<nclus;ic++){
987 AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
988 if(!c)
989 continue;
990 if(!c->IsEMCAL())
991 continue;
5bdfdbea 992 if(c->E()<fMinIsoClusE)
993 continue;
adce44d4 994 Short_t id=-1;
2280356c 995 GetMaxCellEnergy( c, id);
1c86c72c 996 Double_t maxct = cells->GetCellTime(id);
f2a6798b 997 if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc))
1c86c72c 998 continue;
999 Float_t clsPos[3] = {0,0,0};
1000 c->GetPosition(clsPos);
1001 TVector3 cv(clsPos);
9ccace04 1002 cv -= fVecPv;
1c86c72c 1003 Double_t Et = c->E()*TMath::Sin(cv.Theta());
950b5606 1004 Float_t dphi = (cv.Phi()-phicl);
1005 Float_t deta = (cv.Eta()-etacl);
30159e6f 1006 Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1c86c72c 1007 if(R<0.007)
1008 continue;
22814624 1009 if(maxid==id)
1010 continue;
edd083b7 1011 Double_t matchedpt = GetTrackMatchedPt(c->GetTrackMatchedIndex());
950b5606 1012 if(fCpvFromTrack){
1013 if(matchedpt>0 && fRemMatchClus )
1014 continue;
1015 } else {
22603a9b 1016 if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
22603a9b 1017 if(fRemMatchClus){
1018 if(fDebug)
1019 printf("This isolation cluster is matched to a track!++++++++++++++++++++++++++++++++++++++++++++++++++\n");
1020 continue;
1021 }
950b5606 1022 }
1023 }
1c86c72c 1024 Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
1025 if(nEt<0)
1026 printf("nEt=%1.1f\n",nEt);
30159e6f 1027 if(R<fIsoConeR){
1ea1cc08 1028 if(c->GetM02()>0.1 && c->GetM02()<0.3 && !(matchedpt>0)){
1029 TLorentzVector lv, lvec;
1030 lv.SetPtEtaPhiM(Et,cv.Eta(),cv.Phi(),0);
1031 lvec.SetPtEtaPhiM(EtCl,vec.Eta(),vec.Phi(),0);
1032 TLorentzVector lpair = lv + lvec;
e38ddacc 1033 fInConeInvMass += Form("%f;",lpair.M());
1034 if(lpair.M()>0.11 && lpair.M()<0.165){
4ec9dd81 1035 fInConePairedClusEtVsCandEt->Fill(EtCl,Et);
e38ddacc 1036 fInConePairClEt += Form("%f;",Et);
bf10b092 1037 //continue;
e38ddacc 1038 }
1039 else
1040 fInConePairClEt += Form("%f;",0.0);
bf10b092 1041 /*if(lpair.M()>0.52 && lpair.M()<0.58)
1042 continue;*/
1ea1cc08 1043 }
aa119583 1044 totiso += nEt;
30159e6f 1045 if(R<0.04)
1c86c72c 1046 totcore += nEt;
30159e6f 1047 }
1048 else{
1049 if(dphi>fIsoConeR)
1050 continue;
1051 if(deta<fIsoConeR)
1052 continue;
1c86c72c 1053 totband += nEt;
30159e6f 1054 }
1055 }
1056 iso = totiso;
1057 phiband = totband;
1058 core = totcore;
1059}
1060//________________________________________________________________________
1061void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
1062{
bd092f0f 1063 // Get track isolation.
1064
30159e6f 1065 if(!fSelPrimTracks)
1066 return;
f9e2362a 1067 fHigherPtCone = 0;
30159e6f 1068 const Int_t ntracks = fSelPrimTracks->GetEntries();
1069 Double_t totiso=0;
1070 Double_t totband=0;
1071 Double_t totcore=0;
1072 Double_t etacl = vec.Eta();
1073 Double_t phicl = vec.Phi();
1074 if(phicl<0)
1075 phicl+=TMath::TwoPi();
1076 for(int itrack=0;itrack<ntracks;itrack++){
3f4073ba 1077 AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
30159e6f 1078 if(!track)
1079 continue;
1080 Double_t dphi = TMath::Abs(track->Phi()-phicl);
1081 Double_t deta = TMath::Abs(track->Eta()-etacl);
1082 Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1083 Double_t pt = track->Pt();
f9e2362a 1084 if(pt>fHigherPtCone)
1085 fHigherPtCone = pt;
30159e6f 1086 if(R<fIsoConeR){
1087 totiso += track->Pt();
26d2a019 1088 if(R<0.04 && this->fTrCoreRem)
30159e6f 1089 totcore += pt;
1090 }
1091 else{
1092 if(dphi>fIsoConeR)
1093 continue;
1094 if(deta<fIsoConeR)
1095 continue;
1096 totband += track->Pt();
1097 }
1098 }
1099 iso = totiso;
1100 phiband = totband;
1101 core = totcore;
1102}
bd092f0f 1103
30159e6f 1104//________________________________________________________________________
1105Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1106{
1107 // Calculate the energy of cross cells around the leading cell.
950b5606 1108
30159e6f 1109 AliVCaloCells *cells = 0;
2e20efe5 1110 cells = fESDCells;
1111 if (!cells)
1112 cells = fAODCells;
30159e6f 1113 if (!cells)
1114 return 0;
1115
30159e6f 1116 if (!fGeom)
1117 return 0;
1118
1119 Int_t iSupMod = -1;
1120 Int_t iTower = -1;
1121 Int_t iIphi = -1;
1122 Int_t iIeta = -1;
1123 Int_t iphi = -1;
1124 Int_t ieta = -1;
1125 Int_t iphis = -1;
1126 Int_t ietas = -1;
1127
1128 Double_t crossEnergy = 0;
1129
1130 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1131 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1132
1133 Int_t ncells = cluster->GetNCells();
1134 for (Int_t i=0; i<ncells; i++) {
1135 Int_t cellAbsId = cluster->GetCellAbsId(i);
1136 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1137 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1138 Int_t aphidiff = TMath::Abs(iphi-iphis);
1139 if (aphidiff>1)
1140 continue;
1141 Int_t aetadiff = TMath::Abs(ieta-ietas);
1142 if (aetadiff>1)
1143 continue;
1144 if ( (aphidiff==1 && aetadiff==0) ||
1145 (aphidiff==0 && aetadiff==1) ) {
1146 crossEnergy += cells->GetCellAmplitude(cellAbsId);
1147 }
1148 }
1149
1150 return crossEnergy;
1151}
1152
30159e6f 1153//________________________________________________________________________
1154Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1155{
1156 // Get maximum energy of attached cell.
1157
1158 id = -1;
1159
1160 AliVCaloCells *cells = 0;
2e20efe5 1161 cells = fESDCells;
30159e6f 1162 if (!cells)
2e20efe5 1163 cells = fAODCells;
1164 if(!cells)
30159e6f 1165 return 0;
1166
1167 Double_t maxe = 0;
1168 Int_t ncells = cluster->GetNCells();
1169 for (Int_t i=0; i<ncells; i++) {
1170 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1171 if (e>maxe) {
1172 maxe = e;
1173 id = cluster->GetCellAbsId(i);
1174 }
1175 }
1176 return maxe;
1177}
bd092f0f 1178
c7bb0b43 1179//________________________________________________________________________
1180void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
1181{
68faccf4 1182 if(!fStack && !fAODMCParticles)
c7bb0b43 1183 return;
68faccf4 1184 //ESD
1185 if(fStack){
1186 Int_t nTracks = fStack->GetNtrack();
1187 if(fDebug)
1188 printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1189 for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1190 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1191 if(!mcp)
1192 continue;
1193 Int_t pdg = mcp->GetPdgCode();
1194 if(pdg!=22)
1195 continue;
1196 if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1197 continue;
1198 Int_t imom = mcp->GetMother(0);
1199 if(imom<0 || imom>nTracks)
1200 continue;
1201 TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
1202 if(!mcmom)
1203 continue;
1204 Int_t pdgMom = mcmom->GetPdgCode();
1205 Double_t mcphi = mcp->Phi();
1206 Double_t mceta = mcp->Eta();
1207 if((imom==6 || imom==7) && pdgMom==22) {
1208 fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1209 Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1210 fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1211 if(ptsum<2)
1212 fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1213 if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1214 fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1215 if(fNClusForDirPho==0)
1216 fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1217 if(fDebug){
1218 printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
1219 printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
1220 }
1221 }
1222 else{
1223 if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1224 fDecayPhotonPtMC->Fill(mcp->Pt());
40c4f1bf 1225 }
f2acdbe9 1226 }
68faccf4 1227 }
1228 //AOD
1229 else if(fAODMCParticles){
1230 Int_t nTracks = fAODMCParticles->GetEntriesFast();
1231 if(fDebug)
1232 printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1233 for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1234 AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1235 if(!mcp)
1236 continue;
1237 Int_t pdg = mcp->GetPdgCode();
1238 if(pdg!=22)
1239 continue;
1240 if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1241 continue;
1242 Int_t imom = mcp->GetMother();
1243 if(imom<0 || imom>nTracks)
1244 continue;
1245 AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom));
1246 if(!mcmom)
1247 continue;
1248 Int_t pdgMom = mcmom->GetPdgCode();
1249 Double_t mcphi = mcp->Phi();
1250 Double_t mceta = mcp->Eta();
1251 if((imom==6 || imom==7) && pdgMom==22) {
1252 fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1253 Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1254 fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1255 if(ptsum<2)
1256 fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1257 if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1258 fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1259 if(fNClusForDirPho==0)
1260 fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1261 if(fDebug){
1262 printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatus());
1263 printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatus());
1264 }
1265 }
1266 else{
1267 if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1268 fDecayPhotonPtMC->Fill(mcp->Pt());
1269 }
c7bb0b43 1270 }
1271 }
1272}
30159e6f 1273//________________________________________________________________________
f507c09b 1274Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
1275{
1276 if(!c)
84898d7b 1277 return -0.1;
68faccf4 1278 if(!fStack && !fAODMCParticles)
84898d7b 1279 return -0.1;
f507c09b 1280 Int_t clabel = c->GetLabel();
1281 if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
1282 printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
1283 if(!fMcIdFamily.Contains(Form("%d",clabel)))
84898d7b 1284 return -0.1;
f507c09b 1285 fNClusForDirPho++;
1286 TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
1287 Int_t partonpos = partonposstr.Atoi();
1288 if(fDebug)
1289 printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
f507c09b 1290 Float_t clsPos[3] = {0,0,0};
1291 c->GetPosition(clsPos);
1292 TVector3 clsVec(clsPos);
9ccace04 1293 clsVec -= fVecPv;
f507c09b 1294 Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
68faccf4 1295 //ESD
1296 if(fStack){
1297 Int_t nmcp = fStack->GetNtrack();
1298 if(clabel<0 || clabel>nmcp)
1299 return -0.1;
1300 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
1301 if(!mcp)
1302 return -0.1;
1303 if(fDebug){
1304 printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),clabel);
1305 }
1306 Int_t lab1 = mcp->GetFirstDaughter();
1307 if(lab1<0 || lab1>nmcp)
1308 return -0.1;
1309 TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
1310 if(!mcd)
1311 return -0.1;
1312 if(fDebug)
1313 printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->Energy(),mcd->GetPdgCode(),lab1);
1314 if(mcd->GetPdgCode()==22){
1315 fClusEtMcPt->Fill(mcd->Pt(), Et);
1316 fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1317 }
1318 else{
83c10249 1319 if(fDebug)
1320 printf("Warning: daughter of photon parton is not a photon\n");
68faccf4 1321 return -0.1;
1322 }
f507c09b 1323 }
68faccf4 1324 //AOD
1325 else if(fAODMCParticles){
1326 Int_t nmcp = fAODMCParticles->GetEntriesFast();
1327 if(clabel<0 || clabel>nmcp)
1328 return -0.1;
1329 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos));
1330 if(!mcp)
1331 return -0.1;
1332 if(fDebug){
1333 printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),clabel);
1334 }
1335 Int_t lab1 = mcp->GetDaughter(0);
1336 if(lab1<0 || lab1>nmcp)
1337 return -0.1;
1338 AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1));
1339 if(!mcd)
1340 return -0.1;
1341 if(fDebug)
1342 printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->E(),mcd->GetPdgCode(),lab1);
1343 if(mcd->GetPdgCode()==22){
1344 fClusEtMcPt->Fill(mcd->Pt(), Et);
1345 fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1346 }
1347 else{
1348 printf("Warning: daughter of photon parton is not a photon\n");
1349 return -0.1;
1350 }
f507c09b 1351 }
1352 return fDirPhoPt;
1353}
1354//________________________________________________________________________
1355void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
1356{
68faccf4 1357 if(!fStack && !fAODMCParticles)
f507c09b 1358 return;
1359 Int_t selfid = 6;
68faccf4 1360 //ESD
1361 if(fStack){
1362 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
f507c09b 1363 if(!mcp)
1364 return;
68faccf4 1365 if(mcp->GetPdgCode()!=22){
1366 mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
1367 if(!mcp)
1368 return;
1369 }
1370 Int_t daug0f = mcp->GetFirstDaughter();
1371 Int_t daug0l = mcp->GetLastDaughter();
1372 Int_t nd0 = daug0l - daug0f;
1373 if(fDebug)
1374 printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),nd0+1);
1375 fMcIdFamily = Form("%d,",selfid);
1376 GetDaughtersInfo(daug0f,daug0l, selfid,"");
1377 if(fDebug){
1378 printf("\t---- end of the prompt gamma's family tree ----\n\n");
1379 printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1380 }
1381 TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1382 if(!md)
1383 return;
1384 fDirPhoPt = md->Pt();
f507c09b 1385 }
68faccf4 1386 //AOD
1387 else if(fAODMCParticles){
1388 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid));
1389 if(!mcp)
1390 return;
1391 if(mcp->GetPdgCode()!=22){
1392 mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid));
1393 if(!mcp)
1394 return;
1395 }
1396 Int_t daug0f = mcp->GetDaughter(0);
1397 Int_t daug0l = mcp->GetDaughter(mcp->GetNDaughters()-1);
1398 Int_t nd0 = daug0l - daug0f;
1399 if(fDebug)
1400 printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),nd0+1);
1401 fMcIdFamily = Form("%d,",selfid);
1402 GetDaughtersInfo(daug0f,daug0l, selfid,"");
1403 if(fDebug){
1404 printf("\t---- end of the prompt gamma's family tree ----\n\n");
1405 printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1406 }
1407 AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f));
1408 if(!md)
1409 return;
1410 fDirPhoPt = md->Pt();
1411 }
1412
f507c09b 1413}
1414//________________________________________________________________________
22ad7981 1415void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
f507c09b 1416{
68faccf4 1417 if(fStack){
1418 int nmcp = fStack->GetNtrack();
1419 if(firstd<0 || firstd>nmcp)
1420 return;
1421 if(lastd<0 || lastd>nmcp)
1422 return;
1423 if(firstd>lastd){
1424 printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
f507c09b 1425 return;
68faccf4 1426 }
1427 TString indenter = Form("\t%s",inputind);
1428 TParticle *dp = 0x0;
1429 if(fDebug)
1430 printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1431 for(int id=firstd; id<lastd+1; id++){
1432 dp = static_cast<TParticle*>(fStack->Particle(id));
1433 if(!dp)
1434 continue;
1435 Int_t dfd = dp->GetFirstDaughter();
1436 Int_t dld = dp->GetLastDaughter();
1437 Int_t dnd = dld - dfd + 1;
1438 Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1439 if(fDebug)
1440 printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->Energy(), vr, dp->GetPdgCode(), dnd, dfd, dld);
1441 fMcIdFamily += Form("%d,",id);
1442 GetDaughtersInfo(dfd,dld,id,indenter.Data());
1443 }
f507c09b 1444 }
68faccf4 1445 if(fAODMCParticles){
1446 int nmcp = fAODMCParticles->GetEntriesFast();
1447 if(firstd<0 || firstd>nmcp)
1448 return;
1449 if(lastd<0 || lastd>nmcp)
1450 return;
1451 if(firstd>lastd){
1452 printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1453 return;
1454 }
1455 TString indenter = Form("\t%s",inputind);
1456 AliAODMCParticle *dp = 0x0;
f507c09b 1457 if(fDebug)
68faccf4 1458 printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1459 for(int id=firstd; id<lastd+1; id++){
1460 dp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(id));
1461 if(!dp)
1462 continue;
1463 Int_t dfd = dp->GetDaughter(0);
1464 Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1);
1465 Int_t dnd = dld - dfd + 1;
1466 Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv());
1467 if(fDebug)
1468 printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->E(), vr, dp->GetPdgCode(), dnd, dfd, dld);
1469 fMcIdFamily += Form("%d,",id);
1470 GetDaughtersInfo(dfd,dld,id,indenter.Data());
1471 }
f507c09b 1472 }
1473}
f912f9a9 1474
1475//________________________________________________________________________
1476Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1477{
68faccf4 1478 if(!fStack && !fAODMCParticles)
f912f9a9 1479 return 0;
1480 if(fDebug)
1481 printf("Inside GetMcPtSumInCone!!\n");
f912f9a9 1482 Float_t ptsum = 0;
672df96f 1483 TString addpartlabels = "";
68faccf4 1484 //ESD
1485 if(fStack){
1486 Int_t nTracks = fStack->GetNtrack();
1487 for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1488 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1489 if(!mcp)
1490 continue;
1491 Int_t status = mcp->GetStatusCode();
1492 if(status!=1)
1493 continue;
1494 Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1495 if(mcvr>1)
1496 continue;
1497 /*else {
1498 if(fDebug)
f912f9a9 1499 printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
68faccf4 1500 }*/
1501 Float_t dphi = mcp->Phi() - phiclus;
1502 Float_t deta = mcp->Eta() - etaclus;
1503 if(fDebug && TMath::Abs(dphi)<0.01)
1504 printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1505
1506 if(deta>R || dphi>R)
1507 continue;
1508 Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1509 if(dR>R)
1510 continue;
1511 addpartlabels += Form("%d,",iTrack);
1512 if(mcp->Pt()<0.2)
1513 continue;
1514 ptsum += mcp->Pt();
1515 }
1516 }
1517 //AOD
1518 if(fAODMCParticles){
1519 Int_t nTracks = fAODMCParticles->GetEntriesFast();
1520 for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1521 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1522 if(!mcp)
1523 continue;
1524 Int_t status = mcp->GetStatus();
1525 if(status!=1)
1526 continue;
1527 Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv());
1528 if(mcvr>1)
1529 continue;
1530 /*else {
1531 if(fDebug)
1532 printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1533 }*/
1534 Float_t dphi = mcp->Phi() - phiclus;
1535 Float_t deta = mcp->Eta() - etaclus;
1536 if(fDebug && TMath::Abs(dphi)<0.01)
1537 printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1538
1539 if(deta>R || dphi>R)
1540 continue;
1541 Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1542 if(dR>R)
1543 continue;
1544 addpartlabels += Form("%d,",iTrack);
1545 if(mcp->Pt()<0.2)
1546 continue;
1547 ptsum += mcp->Pt();
1548 }
f912f9a9 1549 }
1550 return ptsum;
1551}
f507c09b 1552//________________________________________________________________________
2b7205ad 1553void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1554{
85b52a52 1555
1556 TObjArray *clusters = fESDClusters;
877552c2 1557 //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"),
85b52a52 1558 if (!clusters){
1559 clusters = fAODClusters;
1560 if(fDebug)
1561 printf("ESD clusters empty...");
1562 }
1563 if (!clusters){
1564 if(fDebug)
1565 printf(" and AOD clusters as well!!!\n");
1566 return;
1567 }
2b7205ad 1568 if(!fSelPrimTracks)
1569 return;
1570 const int ntracks = fSelPrimTracks->GetEntriesFast();
112eb594 1571 const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
85b52a52 1572 const Int_t nclus = clusters->GetEntries();
2b7205ad 1573 fNTracks->Fill(ntracks);
1574 fEmcNCells->Fill(ncells);
1575 fEmcNClus->Fill(nclus);
1576 if(fMaxEClus>fECut){
1577 fNTracksECut->Fill(ntracks);
1578 fEmcNCellsCut->Fill(ncells);
1579 fEmcNClusCut->Fill(nclus);
1580 }
1581 for(int it=0;it<ntracks;it++){
1582 AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1583 if(!t)
1584 continue;
e4ba80ad 1585 fTrackPtPhi->Fill(t->Pt(),t->Phi());
1586 fTrackPtEta->Fill(t->Pt(),t->Eta());
1587 if(fMaxEClus>fECut){
1588 fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1589 fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1590 }
ed4488b6 1591 if(t->GetTPCsignal()<80 || t->GetTPCsignal()>100)
1592 continue;
1593 if(t->GetEMCALcluster()<=0 || t->GetEMCALcluster()>nclus)
1594 continue;
1595 AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(t->GetEMCALcluster()));
1596 if(!c)
1597 continue;
1598 fEoverPvsE->Fill(c->E(),c->E()/t->P());
2b7205ad 1599 }
1600 for(int ic=0;ic<nclus;ic++){
85b52a52 1601 AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic));
1602 //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic);
2b7205ad 1603 if(!c)
1604 continue;
3ae97198 1605 if(!c->IsEMCAL())
1606 continue;
e4ba80ad 1607 Float_t clsPos[3] = {0,0,0};
1608 c->GetPosition(clsPos);
1609 TVector3 clsVec(clsPos);
9ccace04 1610 clsVec -= fVecPv;
e4ba80ad 1611 Double_t cphi = clsVec.Phi();
1612 Double_t ceta = clsVec.Eta();
adce44d4 1613 Short_t id = -1;
877552c2 1614 GetMaxCellEnergy( c, id);
134db268 1615 fEmcClusEClusCuts->Fill(c->E(),0);
e4ba80ad 1616 fEmcClusEPhi->Fill(c->E(), cphi);
1617 fEmcClusEEta->Fill(c->E(), ceta);
1618 if(fMaxEClus>fECut){
1619 fEmcClusEPhiCut->Fill(c->E(), cphi);
1620 fEmcClusEEtaCut->Fill(c->E(), ceta);
1621 }
877552c2 1622 Double_t maxt=0;
1623 if(fESDCells)
1624 maxt = fESDCells->GetCellTime(id);
1625 else if(fAODCells)
1626 maxt = fAODCells->GetCellTime(id);
1627 if(IsExotic(c))
1628 continue;
1629 fEmcClusNotExo->Fill(c->E());
134db268 1630 fEmcClusEClusCuts->Fill(c->E(),1);
74d5d8ab 1631 if(fClusIdFromTracks.Contains(Form("%d",ic))){
877552c2 1632 fEmcClusETM2->Fill(c->E());
74d5d8ab 1633 fDetaDphiFromTM->Fill(c->GetTrackDz(),c->GetTrackDx());
1634 }
877552c2 1635 if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1636 fEmcClusETM1->Fill(c->E());
1637 continue;
1638 }
134db268 1639 fEmcClusEClusCuts->Fill(c->E(),2);
992d3f4b 1640 if(TMath::Abs(maxt)>30e-9)
877552c2 1641 continue;
134db268 1642 fEmcClusEClusCuts->Fill(c->E(),3);
877552c2 1643 if(c->GetM02()>0.1)
134db268 1644 fEmcClusEClusCuts->Fill(c->E(),4);
2b7205ad 1645 }
1646}
1647//________________________________________________________________________
1c86c72c 1648Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1649{
1650 Double_t pt = 0;
1651 if(!fTracks)
1652 return pt;
1653 if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1654 if(fDebug)
1655 printf("track-matched index out of track array range!!!\n");
1656 return pt;
1657 }
1658 AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1659 if(!track){
1660 if(fDebug)
1661 printf("track-matched pointer does not exist!!!\n");
1662 return pt;
1663 }
1664 if(fESD){
e1c2239f 1665 if(fPrTrCuts && fPrTrCuts->IsSelected(track))
1666 pt = track->Pt();
1667 else {
1668 if(fCompTrCuts && fCompTrCuts->IsSelected(track))
1669 pt = track->Pt();
1670 else
1671 return pt;
1672 }
1c86c72c 1673 }
1674 return pt;
1675}
1676//________________________________________________________________________
112eb594 1677void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1678{
1679 AliVCaloCells *cells = 0;
1680 cells = fESDCells;
1681 if (!cells)
1682 cells = fAODCells;
1683 if(!cells)
1684 return;
3ae97198 1685 Double_t maxe = 0;
1686 Double_t maxphi = -10;
112eb594 1687 Int_t ncells = cells->GetNumberOfCells();
3ae97198 1688 Double_t eta,phi;
112eb594 1689 for (Int_t i=0; i<ncells; i++) {
1690 Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1691 Double_t e = cells->GetCellAmplitude(absid);
1692 if(e>0.05)
1693 fNCells50++;
3ae97198 1694 else
1695 continue;
1696 fGeom->EtaPhiFromIndex(absid,eta,phi);
1697 if(maxe<e){
1698 maxe = e;
1699 maxphi = phi;
1700 }
112eb594 1701 }
3ae97198 1702 fMaxCellEPhi->Fill(maxe,maxphi);
112eb594 1703
1704}
1705//________________________________________________________________________
34393c41 1706bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c)
1707{
1708 bool isExo = 0;
adce44d4 1709 Short_t id = -1;
34393c41 1710 Double_t Emax = GetMaxCellEnergy( c, id);
1711 Double_t Ecross = GetCrossEnergy( c, id);
1712 if((1-Ecross/Emax)>fExoticCut)
1713 isExo = 1;
1714 return isExo;
1715}
1716//________________________________________________________________________
30159e6f 1717void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1718{
bd092f0f 1719 // Called once at the end of the query.
30159e6f 1720}