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