]> 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
RWGCF converted to native cmake
[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;
7022f12d 836 if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV){
837 TObjArray *inConeInvMassArr = (TObjArray*)fInConeInvMass.Tokenize(";");
e38ddacc 838 TObjArray *inConePairClEt = (TObjArray*)fInConePairClEt.Tokenize(";");
c141c3fd 839 nInConePairs = inConeInvMassArr->GetEntriesFast();
e38ddacc 840 Int_t nInConePi0 = inConePairClEt->GetEntriesFast();
841 if(nInConePairs != nInConePi0)
842 printf("Inconsistent number of in cone pairs!!!\n");
7022f12d 843 for(int ipair=0;ipair<nInConePairs;ipair++){
844 TObjString *obs = (TObjString*)inConeInvMassArr->At(ipair);
e38ddacc 845 TObjString *obet = (TObjString*)inConePairClEt->At(ipair);
7022f12d 846 TString smass = obs->GetString();
e38ddacc 847 TString spairEt = obet->GetString();
7022f12d 848 Double_t pairmass = smass.Atof();
e38ddacc 849 Double_t pairEt = spairEt.Atof();//this must be zero when inv mass outside pi0 range
c141c3fd 850 if(0==ipair && nInConePairs==1)
851 onePairMass = pairmass;
7022f12d 852 if(fDebug)
853 printf("=================+++++++++++++++Inv mass inside the cone for photon range: %1.1f,%1.1f,%1.1f+-++++-+-+-+-++-+-+-\n",Et,pairmass,ceiso+triso);
e38ddacc 854 fEtCandIsoAndIsoWoPairEt->Fill(Et,ceiso+triso,ceiso+triso-pairEt);
7022f12d 855 }
7022f12d 856 }
f9e2362a 857 Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
858 if(Et>10 && Et<15 && dr>0.025){
859 fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
860 if(fDebug)
861 printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
862 }
30159e6f 863 alliso = ceiso + triso;
864 allphiband = cephiband + trphiband;
6f00a63d 865 //allcore = cecore + trcore;
30159e6f 866 Float_t ceisoue = cephiband/phibandArea*netConeArea;
867 Float_t trisoue = trphiband/phibandArea*netConeArea;
868 Float_t allisoue = allphiband/phibandArea*netConeArea;
f912f9a9 869 Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
870 if(fDebug && Et>10)
871 printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
872 if(c->GetM02()>0.5 && c->GetM02()<2.0){
559660d1 873 fMcPtInConeBG->Fill(alliso-allisoue, mcptsum);
874 fMcPtInConeBGnoUE->Fill(alliso, mcptsum);
875 fMcPtInConeTrBGnoUE->Fill(triso, mcptsum);
f912f9a9 876 }
b16ded35 877 if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
559660d1 878 fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum);
879 fMcPtInConeSBGnoUE->Fill(alliso, mcptsum);
880 fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum);
81f2660f 881 if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
559660d1 882 fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue);
883 fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore);
81f2660f 884 }
f912f9a9 885 }
1495ca60 886 if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV)
d834227a 887 fClusEtCPVSBGISO->Fill(Et,alliso - trcore);
1495ca60 888 if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV)
26d2a019 889 fClusEtCPVBGISO->Fill(Et,alliso - trcore);
965c985f 890 const Int_t ndims = fNDimensions;
891 Double_t outputValues[ndims];
e53ab710 892 if(mcptsum<2)
893 ptmc = GetClusSource(c);
894 else
895 ptmc = -0.1;
965c985f 896 outputValues[0] = Et;
897 outputValues[1] = c->GetM02();
1c86c72c 898 outputValues[2] = ceiso/*cecore*/-ceisoue;
965c985f 899 outputValues[3] = triso-trisoue;
1c86c72c 900 outputValues[4] = alliso/*cecore*/-allisoue - trcore;
901 outputValues[5] = ceiso;
902 outputValues[6] = alliso - trcore;
8e574ec0 903 if(fDebug)
904 printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz());
905 if(TMath::Abs(c->GetTrackDx())<0.1)
906 outputValues[7] = c->GetTrackDx();
907 else
908 outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx());
909 if(TMath::Abs(c->GetTrackDz())<0.05)
910 outputValues[8] = c->GetTrackDz();
911 else
912 outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz());
717bb45b 913 outputValues[9] = clsVec.Eta();
914 outputValues[10] = clsVec.Phi();
2e20efe5 915 if(fESDCells)
916 outputValues[11] = fESDCells->GetCellTime(id);
917 else if(fAODCells)
918 outputValues[11] = fAODCells->GetCellTime(id);
16a4050e 919 outputValues[12] = fTrackMult;
f507c09b 920 outputValues[13] = ptmc;
c141c3fd 921 if(nInConePairs == 1)
922 outputValues[14] = onePairMass;
923 else
924 outputValues[14] = -1;
965c985f 925 fHnOutput->Fill(outputValues);
30159e6f 926 if(c->E()>maxE)
927 maxE = c->E();
928 }
929 fNClusHighClusE->Fill(maxE,nclus);
ee3d9c10 930 fMaxEClus = maxE;
0b21f686 931 fNClusEt10->Fill(nclus10);
f507c09b 932 fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
30159e6f 933}
bd092f0f 934
30159e6f 935//________________________________________________________________________
7022f12d 936void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core, Double_t EtCl)
30159e6f 937{
950b5606 938 if(fDebug)
939 printf("....indside GetCeIso funtcion\n");
bd092f0f 940 // Get cell isolation.
2e20efe5 941 AliVCaloCells *cells = fESDCells;
1c86c72c 942 if (!cells){
2e20efe5 943 cells = fAODCells;
1c86c72c 944 if(fDebug)
3aca451b 945 printf("ESD cells empty...\n");
1c86c72c 946 }
947 if (!cells){
948 if(fDebug)
949 printf(" and AOD cells are empty as well!!!\n");
30159e6f 950 return;
1c86c72c 951 }
2e20efe5 952
1c86c72c 953 TObjArray *clusters = fESDClusters;
954 if (!clusters)
955 clusters = fAODClusters;
956 if (!clusters)
957 return;
958
7022f12d 959 fInConeInvMass = "";
e38ddacc 960 fInConePairClEt="";
1c86c72c 961 const Int_t nclus = clusters->GetEntries();
962 //const Int_t ncells = cells->GetNumberOfCells();
30159e6f 963 Float_t totiso=0;
964 Float_t totband=0;
965 Float_t totcore=0;
966 Float_t etacl = vec.Eta();
967 Float_t phicl = vec.Phi();
30159e6f 968 if(phicl<0)
969 phicl+=TMath::TwoPi();
1c86c72c 970 /*Int_t absid = -1;
30159e6f 971 Float_t eta=-1, phi=-1;
972 for(int icell=0;icell<ncells;icell++){
2e20efe5 973 absid = TMath::Abs(cells->GetCellNumber(icell));
974 Float_t celltime = cells->GetCellTime(absid);
f224025f 975 //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
976 if(TMath::Abs(celltime-maxtcl)>2e-8 )
c3f33463 977 continue;
30159e6f 978 if(!fGeom)
979 return;
980 fGeom->EtaPhiFromIndex(absid,eta,phi);
981 Float_t dphi = TMath::Abs(phi-phicl);
982 Float_t deta = TMath::Abs(eta-etacl);
1c86c72c 983 Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
984 for(int ic=0;ic<nclus;ic++){
985 AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
986 if(!c)
987 continue;
988 if(!c->IsEMCAL())
989 continue;
5bdfdbea 990 if(c->E()<fMinIsoClusE)
991 continue;
adce44d4 992 Short_t id=-1;
2280356c 993 GetMaxCellEnergy( c, id);
1c86c72c 994 Double_t maxct = cells->GetCellTime(id);
f2a6798b 995 if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc))
1c86c72c 996 continue;
997 Float_t clsPos[3] = {0,0,0};
998 c->GetPosition(clsPos);
999 TVector3 cv(clsPos);
9ccace04 1000 cv -= fVecPv;
1c86c72c 1001 Double_t Et = c->E()*TMath::Sin(cv.Theta());
950b5606 1002 Float_t dphi = (cv.Phi()-phicl);
1003 Float_t deta = (cv.Eta()-etacl);
30159e6f 1004 Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1c86c72c 1005 if(R<0.007)
1006 continue;
22814624 1007 if(maxid==id)
1008 continue;
edd083b7 1009 Double_t matchedpt = GetTrackMatchedPt(c->GetTrackMatchedIndex());
950b5606 1010 if(fCpvFromTrack){
1011 if(matchedpt>0 && fRemMatchClus )
1012 continue;
1013 } else {
22603a9b 1014 if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
22603a9b 1015 if(fRemMatchClus){
1016 if(fDebug)
1017 printf("This isolation cluster is matched to a track!++++++++++++++++++++++++++++++++++++++++++++++++++\n");
1018 continue;
1019 }
950b5606 1020 }
1021 }
1c86c72c 1022 Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
1023 if(nEt<0)
1024 printf("nEt=%1.1f\n",nEt);
30159e6f 1025 if(R<fIsoConeR){
1ea1cc08 1026 if(c->GetM02()>0.1 && c->GetM02()<0.3 && !(matchedpt>0)){
1027 TLorentzVector lv, lvec;
1028 lv.SetPtEtaPhiM(Et,cv.Eta(),cv.Phi(),0);
1029 lvec.SetPtEtaPhiM(EtCl,vec.Eta(),vec.Phi(),0);
1030 TLorentzVector lpair = lv + lvec;
e38ddacc 1031 fInConeInvMass += Form("%f;",lpair.M());
1032 if(lpair.M()>0.11 && lpair.M()<0.165){
4ec9dd81 1033 fInConePairedClusEtVsCandEt->Fill(EtCl,Et);
e38ddacc 1034 fInConePairClEt += Form("%f;",Et);
aa119583 1035 continue;
e38ddacc 1036 }
1037 else
1038 fInConePairClEt += Form("%f;",0.0);
aa119583 1039 if(lpair.M()>0.52 && lpair.M()<0.58)
1040 continue;
1ea1cc08 1041 }
aa119583 1042 totiso += nEt;
30159e6f 1043 if(R<0.04)
1c86c72c 1044 totcore += nEt;
30159e6f 1045 }
1046 else{
1047 if(dphi>fIsoConeR)
1048 continue;
1049 if(deta<fIsoConeR)
1050 continue;
1c86c72c 1051 totband += nEt;
30159e6f 1052 }
1053 }
1054 iso = totiso;
1055 phiband = totband;
1056 core = totcore;
1057}
1058//________________________________________________________________________
1059void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
1060{
bd092f0f 1061 // Get track isolation.
1062
30159e6f 1063 if(!fSelPrimTracks)
1064 return;
f9e2362a 1065 fHigherPtCone = 0;
30159e6f 1066 const Int_t ntracks = fSelPrimTracks->GetEntries();
1067 Double_t totiso=0;
1068 Double_t totband=0;
1069 Double_t totcore=0;
1070 Double_t etacl = vec.Eta();
1071 Double_t phicl = vec.Phi();
1072 if(phicl<0)
1073 phicl+=TMath::TwoPi();
1074 for(int itrack=0;itrack<ntracks;itrack++){
3f4073ba 1075 AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
30159e6f 1076 if(!track)
1077 continue;
1078 Double_t dphi = TMath::Abs(track->Phi()-phicl);
1079 Double_t deta = TMath::Abs(track->Eta()-etacl);
1080 Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1081 Double_t pt = track->Pt();
f9e2362a 1082 if(pt>fHigherPtCone)
1083 fHigherPtCone = pt;
30159e6f 1084 if(R<fIsoConeR){
1085 totiso += track->Pt();
26d2a019 1086 if(R<0.04 && this->fTrCoreRem)
30159e6f 1087 totcore += pt;
1088 }
1089 else{
1090 if(dphi>fIsoConeR)
1091 continue;
1092 if(deta<fIsoConeR)
1093 continue;
1094 totband += track->Pt();
1095 }
1096 }
1097 iso = totiso;
1098 phiband = totband;
1099 core = totcore;
1100}
bd092f0f 1101
30159e6f 1102//________________________________________________________________________
1103Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1104{
1105 // Calculate the energy of cross cells around the leading cell.
950b5606 1106
30159e6f 1107 AliVCaloCells *cells = 0;
2e20efe5 1108 cells = fESDCells;
1109 if (!cells)
1110 cells = fAODCells;
30159e6f 1111 if (!cells)
1112 return 0;
1113
30159e6f 1114 if (!fGeom)
1115 return 0;
1116
1117 Int_t iSupMod = -1;
1118 Int_t iTower = -1;
1119 Int_t iIphi = -1;
1120 Int_t iIeta = -1;
1121 Int_t iphi = -1;
1122 Int_t ieta = -1;
1123 Int_t iphis = -1;
1124 Int_t ietas = -1;
1125
1126 Double_t crossEnergy = 0;
1127
1128 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1129 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1130
1131 Int_t ncells = cluster->GetNCells();
1132 for (Int_t i=0; i<ncells; i++) {
1133 Int_t cellAbsId = cluster->GetCellAbsId(i);
1134 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1135 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1136 Int_t aphidiff = TMath::Abs(iphi-iphis);
1137 if (aphidiff>1)
1138 continue;
1139 Int_t aetadiff = TMath::Abs(ieta-ietas);
1140 if (aetadiff>1)
1141 continue;
1142 if ( (aphidiff==1 && aetadiff==0) ||
1143 (aphidiff==0 && aetadiff==1) ) {
1144 crossEnergy += cells->GetCellAmplitude(cellAbsId);
1145 }
1146 }
1147
1148 return crossEnergy;
1149}
1150
30159e6f 1151//________________________________________________________________________
1152Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1153{
1154 // Get maximum energy of attached cell.
1155
1156 id = -1;
1157
1158 AliVCaloCells *cells = 0;
2e20efe5 1159 cells = fESDCells;
30159e6f 1160 if (!cells)
2e20efe5 1161 cells = fAODCells;
1162 if(!cells)
30159e6f 1163 return 0;
1164
1165 Double_t maxe = 0;
1166 Int_t ncells = cluster->GetNCells();
1167 for (Int_t i=0; i<ncells; i++) {
1168 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1169 if (e>maxe) {
1170 maxe = e;
1171 id = cluster->GetCellAbsId(i);
1172 }
1173 }
1174 return maxe;
1175}
bd092f0f 1176
c7bb0b43 1177//________________________________________________________________________
1178void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
1179{
68faccf4 1180 if(!fStack && !fAODMCParticles)
c7bb0b43 1181 return;
68faccf4 1182 //ESD
1183 if(fStack){
1184 Int_t nTracks = fStack->GetNtrack();
1185 if(fDebug)
1186 printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1187 for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1188 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1189 if(!mcp)
1190 continue;
1191 Int_t pdg = mcp->GetPdgCode();
1192 if(pdg!=22)
1193 continue;
1194 if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1195 continue;
1196 Int_t imom = mcp->GetMother(0);
1197 if(imom<0 || imom>nTracks)
1198 continue;
1199 TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
1200 if(!mcmom)
1201 continue;
1202 Int_t pdgMom = mcmom->GetPdgCode();
1203 Double_t mcphi = mcp->Phi();
1204 Double_t mceta = mcp->Eta();
1205 if((imom==6 || imom==7) && pdgMom==22) {
1206 fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1207 Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1208 fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1209 if(ptsum<2)
1210 fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1211 if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1212 fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1213 if(fNClusForDirPho==0)
1214 fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1215 if(fDebug){
1216 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());
1217 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());
1218 }
1219 }
1220 else{
1221 if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1222 fDecayPhotonPtMC->Fill(mcp->Pt());
40c4f1bf 1223 }
f2acdbe9 1224 }
68faccf4 1225 }
1226 //AOD
1227 else if(fAODMCParticles){
1228 Int_t nTracks = fAODMCParticles->GetEntriesFast();
1229 if(fDebug)
1230 printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1231 for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1232 AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1233 if(!mcp)
1234 continue;
1235 Int_t pdg = mcp->GetPdgCode();
1236 if(pdg!=22)
1237 continue;
1238 if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1239 continue;
1240 Int_t imom = mcp->GetMother();
1241 if(imom<0 || imom>nTracks)
1242 continue;
1243 AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom));
1244 if(!mcmom)
1245 continue;
1246 Int_t pdgMom = mcmom->GetPdgCode();
1247 Double_t mcphi = mcp->Phi();
1248 Double_t mceta = mcp->Eta();
1249 if((imom==6 || imom==7) && pdgMom==22) {
1250 fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1251 Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1252 fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1253 if(ptsum<2)
1254 fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1255 if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1256 fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1257 if(fNClusForDirPho==0)
1258 fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1259 if(fDebug){
1260 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());
1261 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());
1262 }
1263 }
1264 else{
1265 if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1266 fDecayPhotonPtMC->Fill(mcp->Pt());
1267 }
c7bb0b43 1268 }
1269 }
1270}
30159e6f 1271//________________________________________________________________________
f507c09b 1272Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
1273{
1274 if(!c)
84898d7b 1275 return -0.1;
68faccf4 1276 if(!fStack && !fAODMCParticles)
84898d7b 1277 return -0.1;
f507c09b 1278 Int_t clabel = c->GetLabel();
1279 if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
1280 printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
1281 if(!fMcIdFamily.Contains(Form("%d",clabel)))
84898d7b 1282 return -0.1;
f507c09b 1283 fNClusForDirPho++;
1284 TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
1285 Int_t partonpos = partonposstr.Atoi();
1286 if(fDebug)
1287 printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
f507c09b 1288 Float_t clsPos[3] = {0,0,0};
1289 c->GetPosition(clsPos);
1290 TVector3 clsVec(clsPos);
9ccace04 1291 clsVec -= fVecPv;
f507c09b 1292 Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
68faccf4 1293 //ESD
1294 if(fStack){
1295 Int_t nmcp = fStack->GetNtrack();
1296 if(clabel<0 || clabel>nmcp)
1297 return -0.1;
1298 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
1299 if(!mcp)
1300 return -0.1;
1301 if(fDebug){
1302 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);
1303 }
1304 Int_t lab1 = mcp->GetFirstDaughter();
1305 if(lab1<0 || lab1>nmcp)
1306 return -0.1;
1307 TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
1308 if(!mcd)
1309 return -0.1;
1310 if(fDebug)
1311 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);
1312 if(mcd->GetPdgCode()==22){
1313 fClusEtMcPt->Fill(mcd->Pt(), Et);
1314 fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1315 }
1316 else{
83c10249 1317 if(fDebug)
1318 printf("Warning: daughter of photon parton is not a photon\n");
68faccf4 1319 return -0.1;
1320 }
f507c09b 1321 }
68faccf4 1322 //AOD
1323 else if(fAODMCParticles){
1324 Int_t nmcp = fAODMCParticles->GetEntriesFast();
1325 if(clabel<0 || clabel>nmcp)
1326 return -0.1;
1327 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos));
1328 if(!mcp)
1329 return -0.1;
1330 if(fDebug){
1331 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);
1332 }
1333 Int_t lab1 = mcp->GetDaughter(0);
1334 if(lab1<0 || lab1>nmcp)
1335 return -0.1;
1336 AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1));
1337 if(!mcd)
1338 return -0.1;
1339 if(fDebug)
1340 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);
1341 if(mcd->GetPdgCode()==22){
1342 fClusEtMcPt->Fill(mcd->Pt(), Et);
1343 fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1344 }
1345 else{
1346 printf("Warning: daughter of photon parton is not a photon\n");
1347 return -0.1;
1348 }
f507c09b 1349 }
1350 return fDirPhoPt;
1351}
1352//________________________________________________________________________
1353void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
1354{
68faccf4 1355 if(!fStack && !fAODMCParticles)
f507c09b 1356 return;
1357 Int_t selfid = 6;
68faccf4 1358 //ESD
1359 if(fStack){
1360 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
f507c09b 1361 if(!mcp)
1362 return;
68faccf4 1363 if(mcp->GetPdgCode()!=22){
1364 mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
1365 if(!mcp)
1366 return;
1367 }
1368 Int_t daug0f = mcp->GetFirstDaughter();
1369 Int_t daug0l = mcp->GetLastDaughter();
1370 Int_t nd0 = daug0l - daug0f;
1371 if(fDebug)
1372 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);
1373 fMcIdFamily = Form("%d,",selfid);
1374 GetDaughtersInfo(daug0f,daug0l, selfid,"");
1375 if(fDebug){
1376 printf("\t---- end of the prompt gamma's family tree ----\n\n");
1377 printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1378 }
1379 TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1380 if(!md)
1381 return;
1382 fDirPhoPt = md->Pt();
f507c09b 1383 }
68faccf4 1384 //AOD
1385 else if(fAODMCParticles){
1386 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid));
1387 if(!mcp)
1388 return;
1389 if(mcp->GetPdgCode()!=22){
1390 mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid));
1391 if(!mcp)
1392 return;
1393 }
1394 Int_t daug0f = mcp->GetDaughter(0);
1395 Int_t daug0l = mcp->GetDaughter(mcp->GetNDaughters()-1);
1396 Int_t nd0 = daug0l - daug0f;
1397 if(fDebug)
1398 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);
1399 fMcIdFamily = Form("%d,",selfid);
1400 GetDaughtersInfo(daug0f,daug0l, selfid,"");
1401 if(fDebug){
1402 printf("\t---- end of the prompt gamma's family tree ----\n\n");
1403 printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1404 }
1405 AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f));
1406 if(!md)
1407 return;
1408 fDirPhoPt = md->Pt();
1409 }
1410
f507c09b 1411}
1412//________________________________________________________________________
22ad7981 1413void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
f507c09b 1414{
68faccf4 1415 if(fStack){
1416 int nmcp = fStack->GetNtrack();
1417 if(firstd<0 || firstd>nmcp)
1418 return;
1419 if(lastd<0 || lastd>nmcp)
1420 return;
1421 if(firstd>lastd){
1422 printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
f507c09b 1423 return;
68faccf4 1424 }
1425 TString indenter = Form("\t%s",inputind);
1426 TParticle *dp = 0x0;
1427 if(fDebug)
1428 printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1429 for(int id=firstd; id<lastd+1; id++){
1430 dp = static_cast<TParticle*>(fStack->Particle(id));
1431 if(!dp)
1432 continue;
1433 Int_t dfd = dp->GetFirstDaughter();
1434 Int_t dld = dp->GetLastDaughter();
1435 Int_t dnd = dld - dfd + 1;
1436 Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1437 if(fDebug)
1438 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);
1439 fMcIdFamily += Form("%d,",id);
1440 GetDaughtersInfo(dfd,dld,id,indenter.Data());
1441 }
f507c09b 1442 }
68faccf4 1443 if(fAODMCParticles){
1444 int nmcp = fAODMCParticles->GetEntriesFast();
1445 if(firstd<0 || firstd>nmcp)
1446 return;
1447 if(lastd<0 || lastd>nmcp)
1448 return;
1449 if(firstd>lastd){
1450 printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1451 return;
1452 }
1453 TString indenter = Form("\t%s",inputind);
1454 AliAODMCParticle *dp = 0x0;
f507c09b 1455 if(fDebug)
68faccf4 1456 printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1457 for(int id=firstd; id<lastd+1; id++){
1458 dp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(id));
1459 if(!dp)
1460 continue;
1461 Int_t dfd = dp->GetDaughter(0);
1462 Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1);
1463 Int_t dnd = dld - dfd + 1;
1464 Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv());
1465 if(fDebug)
1466 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);
1467 fMcIdFamily += Form("%d,",id);
1468 GetDaughtersInfo(dfd,dld,id,indenter.Data());
1469 }
f507c09b 1470 }
1471}
f912f9a9 1472
1473//________________________________________________________________________
1474Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1475{
68faccf4 1476 if(!fStack && !fAODMCParticles)
f912f9a9 1477 return 0;
1478 if(fDebug)
1479 printf("Inside GetMcPtSumInCone!!\n");
f912f9a9 1480 Float_t ptsum = 0;
672df96f 1481 TString addpartlabels = "";
68faccf4 1482 //ESD
1483 if(fStack){
1484 Int_t nTracks = fStack->GetNtrack();
1485 for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1486 TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1487 if(!mcp)
1488 continue;
1489 Int_t status = mcp->GetStatusCode();
1490 if(status!=1)
1491 continue;
1492 Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1493 if(mcvr>1)
1494 continue;
1495 /*else {
1496 if(fDebug)
f912f9a9 1497 printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
68faccf4 1498 }*/
1499 Float_t dphi = mcp->Phi() - phiclus;
1500 Float_t deta = mcp->Eta() - etaclus;
1501 if(fDebug && TMath::Abs(dphi)<0.01)
1502 printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1503
1504 if(deta>R || dphi>R)
1505 continue;
1506 Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1507 if(dR>R)
1508 continue;
1509 addpartlabels += Form("%d,",iTrack);
1510 if(mcp->Pt()<0.2)
1511 continue;
1512 ptsum += mcp->Pt();
1513 }
1514 }
1515 //AOD
1516 if(fAODMCParticles){
1517 Int_t nTracks = fAODMCParticles->GetEntriesFast();
1518 for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1519 AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1520 if(!mcp)
1521 continue;
1522 Int_t status = mcp->GetStatus();
1523 if(status!=1)
1524 continue;
1525 Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv());
1526 if(mcvr>1)
1527 continue;
1528 /*else {
1529 if(fDebug)
1530 printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1531 }*/
1532 Float_t dphi = mcp->Phi() - phiclus;
1533 Float_t deta = mcp->Eta() - etaclus;
1534 if(fDebug && TMath::Abs(dphi)<0.01)
1535 printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1536
1537 if(deta>R || dphi>R)
1538 continue;
1539 Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1540 if(dR>R)
1541 continue;
1542 addpartlabels += Form("%d,",iTrack);
1543 if(mcp->Pt()<0.2)
1544 continue;
1545 ptsum += mcp->Pt();
1546 }
f912f9a9 1547 }
1548 return ptsum;
1549}
f507c09b 1550//________________________________________________________________________
2b7205ad 1551void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1552{
85b52a52 1553
1554 TObjArray *clusters = fESDClusters;
877552c2 1555 //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"),
85b52a52 1556 if (!clusters){
1557 clusters = fAODClusters;
1558 if(fDebug)
1559 printf("ESD clusters empty...");
1560 }
1561 if (!clusters){
1562 if(fDebug)
1563 printf(" and AOD clusters as well!!!\n");
1564 return;
1565 }
2b7205ad 1566 if(!fSelPrimTracks)
1567 return;
1568 const int ntracks = fSelPrimTracks->GetEntriesFast();
112eb594 1569 const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
85b52a52 1570 const Int_t nclus = clusters->GetEntries();
2b7205ad 1571 fNTracks->Fill(ntracks);
1572 fEmcNCells->Fill(ncells);
1573 fEmcNClus->Fill(nclus);
1574 if(fMaxEClus>fECut){
1575 fNTracksECut->Fill(ntracks);
1576 fEmcNCellsCut->Fill(ncells);
1577 fEmcNClusCut->Fill(nclus);
1578 }
1579 for(int it=0;it<ntracks;it++){
1580 AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1581 if(!t)
1582 continue;
e4ba80ad 1583 fTrackPtPhi->Fill(t->Pt(),t->Phi());
1584 fTrackPtEta->Fill(t->Pt(),t->Eta());
1585 if(fMaxEClus>fECut){
1586 fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1587 fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1588 }
ed4488b6 1589 if(t->GetTPCsignal()<80 || t->GetTPCsignal()>100)
1590 continue;
1591 if(t->GetEMCALcluster()<=0 || t->GetEMCALcluster()>nclus)
1592 continue;
1593 AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(t->GetEMCALcluster()));
1594 if(!c)
1595 continue;
1596 fEoverPvsE->Fill(c->E(),c->E()/t->P());
2b7205ad 1597 }
1598 for(int ic=0;ic<nclus;ic++){
85b52a52 1599 AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic));
1600 //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic);
2b7205ad 1601 if(!c)
1602 continue;
3ae97198 1603 if(!c->IsEMCAL())
1604 continue;
e4ba80ad 1605 Float_t clsPos[3] = {0,0,0};
1606 c->GetPosition(clsPos);
1607 TVector3 clsVec(clsPos);
9ccace04 1608 clsVec -= fVecPv;
e4ba80ad 1609 Double_t cphi = clsVec.Phi();
1610 Double_t ceta = clsVec.Eta();
adce44d4 1611 Short_t id = -1;
877552c2 1612 GetMaxCellEnergy( c, id);
134db268 1613 fEmcClusEClusCuts->Fill(c->E(),0);
e4ba80ad 1614 fEmcClusEPhi->Fill(c->E(), cphi);
1615 fEmcClusEEta->Fill(c->E(), ceta);
1616 if(fMaxEClus>fECut){
1617 fEmcClusEPhiCut->Fill(c->E(), cphi);
1618 fEmcClusEEtaCut->Fill(c->E(), ceta);
1619 }
877552c2 1620 Double_t maxt=0;
1621 if(fESDCells)
1622 maxt = fESDCells->GetCellTime(id);
1623 else if(fAODCells)
1624 maxt = fAODCells->GetCellTime(id);
1625 if(IsExotic(c))
1626 continue;
1627 fEmcClusNotExo->Fill(c->E());
134db268 1628 fEmcClusEClusCuts->Fill(c->E(),1);
74d5d8ab 1629 if(fClusIdFromTracks.Contains(Form("%d",ic))){
877552c2 1630 fEmcClusETM2->Fill(c->E());
74d5d8ab 1631 fDetaDphiFromTM->Fill(c->GetTrackDz(),c->GetTrackDx());
1632 }
877552c2 1633 if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1634 fEmcClusETM1->Fill(c->E());
1635 continue;
1636 }
134db268 1637 fEmcClusEClusCuts->Fill(c->E(),2);
992d3f4b 1638 if(TMath::Abs(maxt)>30e-9)
877552c2 1639 continue;
134db268 1640 fEmcClusEClusCuts->Fill(c->E(),3);
877552c2 1641 if(c->GetM02()>0.1)
134db268 1642 fEmcClusEClusCuts->Fill(c->E(),4);
2b7205ad 1643 }
1644}
1645//________________________________________________________________________
1c86c72c 1646Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1647{
1648 Double_t pt = 0;
1649 if(!fTracks)
1650 return pt;
1651 if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1652 if(fDebug)
1653 printf("track-matched index out of track array range!!!\n");
1654 return pt;
1655 }
1656 AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1657 if(!track){
1658 if(fDebug)
1659 printf("track-matched pointer does not exist!!!\n");
1660 return pt;
1661 }
1662 if(fESD){
c45212c6 1663 if(!fPrTrCuts && !fCompTrCuts)
1c86c72c 1664 return pt;
c45212c6 1665 if(!fPrTrCuts->IsSelected(track) && !fCompTrCuts->IsSelected(track))
1c86c72c 1666 return pt;
1667 pt = track->Pt();
1668 }
1669 return pt;
1670}
1671//________________________________________________________________________
112eb594 1672void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1673{
1674 AliVCaloCells *cells = 0;
1675 cells = fESDCells;
1676 if (!cells)
1677 cells = fAODCells;
1678 if(!cells)
1679 return;
3ae97198 1680 Double_t maxe = 0;
1681 Double_t maxphi = -10;
112eb594 1682 Int_t ncells = cells->GetNumberOfCells();
3ae97198 1683 Double_t eta,phi;
112eb594 1684 for (Int_t i=0; i<ncells; i++) {
1685 Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1686 Double_t e = cells->GetCellAmplitude(absid);
1687 if(e>0.05)
1688 fNCells50++;
3ae97198 1689 else
1690 continue;
1691 fGeom->EtaPhiFromIndex(absid,eta,phi);
1692 if(maxe<e){
1693 maxe = e;
1694 maxphi = phi;
1695 }
112eb594 1696 }
3ae97198 1697 fMaxCellEPhi->Fill(maxe,maxphi);
112eb594 1698
1699}
1700//________________________________________________________________________
34393c41 1701bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c)
1702{
1703 bool isExo = 0;
adce44d4 1704 Short_t id = -1;
34393c41 1705 Double_t Emax = GetMaxCellEnergy( c, id);
1706 Double_t Ecross = GetCrossEnergy( c, id);
1707 if((1-Ecross/Emax)>fExoticCut)
1708 isExo = 1;
1709 return isExo;
1710}
1711//________________________________________________________________________
30159e6f 1712void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1713{
bd092f0f 1714 // Called once at the end of the query.
30159e6f 1715}