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