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