]>
Commit | Line | Data |
---|---|---|
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" | |
2e20efe5 | 32 | #include "AliVEvent.h" |
33 | #include "AliVTrack.h" | |
30159e6f | 34 | #include "AliV0vertexer.h" |
30159e6f | 35 | #include "AliVCluster.h" |
822c9944 | 36 | #include "AliOADBContainer.h" |
30159e6f | 37 | |
c2aad3ae | 38 | #include <iostream> |
39 | using std::cout; | |
40 | using std::endl; | |
41 | ||
30159e6f | 42 | ClassImp(AliAnalysisTaskEMCALIsoPhoton) |
43 | ||
44 | //________________________________________________________________________ | |
bd092f0f | 45 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : |
46 | AliAnalysisTaskSE(), | |
2e20efe5 | 47 | fESDClusters(0), |
48 | fAODClusters(0), | |
bd092f0f | 49 | fSelPrimTracks(0), |
3f4073ba | 50 | fTracks(0), |
2e20efe5 | 51 | fESDCells(0), |
52 | fAODCells(0), | |
bd092f0f | 53 | fPrTrCuts(0), |
54 | fGeom(0x0), | |
9148fa89 | 55 | fGeoName("EMCAL_COMPLETEV1"), |
4d4ce2d2 | 56 | fOADBContainer(0), |
9148fa89 | 57 | fPeriod("LHC11c"), |
58 | fTrigBit("kEMC7"), | |
59 | fIsTrain(0), | |
c7bb0b43 | 60 | fIsMc(0), |
ecd47673 | 61 | fDebug(0), |
f3843637 | 62 | fPathStrOpt("/"), |
9148fa89 | 63 | fExoticCut(0.97), |
64 | fIsoConeR(0.4), | |
965c985f | 65 | fNDimensions(7), |
66 | fECut(3.), | |
16a4050e | 67 | fTrackMult(0), |
f507c09b | 68 | fMcIdFamily(""), |
69 | fNClusForDirPho(0), | |
70 | fDirPhoPt(0), | |
f9e2362a | 71 | fHigherPtCone(0), |
26a8db0f | 72 | fImportGeometryFromFile(0), |
73 | fImportGeometryFilePath(""), | |
2b7205ad | 74 | fMaxPtTrack(0), |
75 | fMaxEClus(0), | |
112eb594 | 76 | fNCells50(0), |
48aab590 | 77 | fFilterBit(0), |
bd092f0f | 78 | fESD(0), |
2e20efe5 | 79 | fAOD(0), |
c7bb0b43 | 80 | fMCEvent(0), |
81 | fStack(0), | |
bd092f0f | 82 | fOutputList(0), |
83 | fEvtSel(0), | |
0b21f686 | 84 | fNClusEt10(0), |
cc57d293 | 85 | fRecoPV(0), |
f507c09b | 86 | fPVtxZ(0), |
87 | fTrMultDist(0), | |
caaf99d3 | 88 | fMCDirPhotonPtEtaPhi(0), |
e53ab710 | 89 | fMCIsoDirPhotonPtEtaPhi(0), |
c7bb0b43 | 90 | fDecayPhotonPtMC(0), |
bd092f0f | 91 | fCellAbsIdVsAmpl(0), |
16a4050e | 92 | fNClusHighClusE(0), |
f9e2362a | 93 | fHigherPtConeM02(0), |
f507c09b | 94 | fClusEtMcPt(0), |
95 | fClusMcDetaDphi(0), | |
96 | fNClusPerPho(0), | |
f912f9a9 | 97 | fMcPtInConeBG(0), |
98 | fMcPtInConeSBG(0), | |
99 | fMcPtInConeBGnoUE(0), | |
100 | fMcPtInConeSBGnoUE(0), | |
81f2660f | 101 | fAllIsoEtMcGamma(0), |
102 | fAllIsoNoUeEtMcGamma(0), | |
092ceec8 | 103 | fMCDirPhotonPtEtaPhiNoClus(0), |
2b7205ad | 104 | fHnOutput(0), |
105 | fQAList(0), | |
106 | fNTracks(0), | |
107 | fEmcNCells(0), | |
108 | fEmcNClus(0), | |
109 | fEmcNClusCut(0), | |
110 | fNTracksECut(0), | |
111 | fEmcNCellsCut(0), | |
e4ba80ad | 112 | fEmcClusEPhi(0), |
113 | fEmcClusEPhiCut(0), | |
114 | fEmcClusEEta(0), | |
115 | fEmcClusEEtaCut(0), | |
116 | fTrackPtPhi(0), | |
117 | fTrackPtPhiCut(0), | |
118 | fTrackPtEta(0), | |
3ae97198 | 119 | fTrackPtEtaCut(0), |
120 | fMaxCellEPhi(0) | |
bd092f0f | 121 | { |
122 | // Default constructor. | |
26a8db0f | 123 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
bd092f0f | 124 | } |
125 | ||
126 | //________________________________________________________________________ | |
127 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : | |
128 | AliAnalysisTaskSE(name), | |
2e20efe5 | 129 | fESDClusters(0), |
130 | fAODClusters(0), | |
30159e6f | 131 | fSelPrimTracks(0), |
3f4073ba | 132 | fTracks(0), |
2e20efe5 | 133 | fESDCells(0), |
134 | fAODCells(0), | |
30159e6f | 135 | fPrTrCuts(0), |
136 | fGeom(0x0), | |
137 | fGeoName("EMCAL_COMPLETEV1"), | |
4d4ce2d2 | 138 | fOADBContainer(0), |
30159e6f | 139 | fPeriod("LHC11c"), |
751194e8 | 140 | fTrigBit("kEMC7"), |
30159e6f | 141 | fIsTrain(0), |
c7bb0b43 | 142 | fIsMc(0), |
ecd47673 | 143 | fDebug(0), |
f3843637 | 144 | fPathStrOpt("/"), |
30159e6f | 145 | fExoticCut(0.97), |
146 | fIsoConeR(0.4), | |
965c985f | 147 | fNDimensions(7), |
148 | fECut(3.), | |
16a4050e | 149 | fTrackMult(0), |
f507c09b | 150 | fMcIdFamily(""), |
151 | fNClusForDirPho(0), | |
152 | fDirPhoPt(0), | |
f9e2362a | 153 | fHigherPtCone(0), |
26a8db0f | 154 | fImportGeometryFromFile(0), |
155 | fImportGeometryFilePath(""), | |
2b7205ad | 156 | fMaxPtTrack(0), |
157 | fMaxEClus(0), | |
112eb594 | 158 | fNCells50(0), |
48aab590 | 159 | fFilterBit(0), |
30159e6f | 160 | fESD(0), |
2e20efe5 | 161 | fAOD(0), |
c7bb0b43 | 162 | fMCEvent(0), |
163 | fStack(0), | |
30159e6f | 164 | fOutputList(0), |
30159e6f | 165 | fEvtSel(0), |
0b21f686 | 166 | fNClusEt10(0), |
cc57d293 | 167 | fRecoPV(0), |
bd092f0f | 168 | fPVtxZ(0), |
f507c09b | 169 | fTrMultDist(0), |
caaf99d3 | 170 | fMCDirPhotonPtEtaPhi(0), |
e53ab710 | 171 | fMCIsoDirPhotonPtEtaPhi(0), |
c7bb0b43 | 172 | fDecayPhotonPtMC(0), |
173 | fCellAbsIdVsAmpl(0), | |
bd092f0f | 174 | fNClusHighClusE(0), |
f9e2362a | 175 | fHigherPtConeM02(0), |
f507c09b | 176 | fClusEtMcPt(0), |
177 | fClusMcDetaDphi(0), | |
178 | fNClusPerPho(0), | |
f912f9a9 | 179 | fMcPtInConeBG(0), |
180 | fMcPtInConeSBG(0), | |
181 | fMcPtInConeBGnoUE(0), | |
182 | fMcPtInConeSBGnoUE(0), | |
81f2660f | 183 | fAllIsoEtMcGamma(0), |
184 | fAllIsoNoUeEtMcGamma(0), | |
092ceec8 | 185 | fMCDirPhotonPtEtaPhiNoClus(0), |
2b7205ad | 186 | fHnOutput(0), |
187 | fQAList(0), | |
188 | fNTracks(0), | |
189 | fEmcNCells(0), | |
190 | fEmcNClus(0), | |
191 | fEmcNClusCut(0), | |
192 | fNTracksECut(0), | |
193 | fEmcNCellsCut(0), | |
e4ba80ad | 194 | fEmcClusEPhi(0), |
195 | fEmcClusEPhiCut(0), | |
196 | fEmcClusEEta(0), | |
197 | fEmcClusEEtaCut(0), | |
198 | fTrackPtPhi(0), | |
199 | fTrackPtPhiCut(0), | |
200 | fTrackPtEta(0), | |
3ae97198 | 201 | fTrackPtEtaCut(0), |
202 | fMaxCellEPhi(0) | |
30159e6f | 203 | { |
204 | // Constructor | |
205 | ||
206 | // Define input and output slots here | |
207 | // Input slot #0 works with a TChain | |
208 | DefineInput(0, TChain::Class()); | |
209 | // Output slot #0 id reserved by the base class for AOD | |
210 | // Output slot #1 writes into a TH1 container | |
211 | DefineOutput(1, TList::Class()); | |
2b7205ad | 212 | DefineOutput(2, TList::Class()); |
30159e6f | 213 | } |
bd092f0f | 214 | |
30159e6f | 215 | //________________________________________________________________________ |
216 | void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() | |
217 | { | |
bd092f0f | 218 | // Create histograms, called once. |
30159e6f | 219 | |
2e20efe5 | 220 | fESDClusters = new TObjArray(); |
30159e6f | 221 | fSelPrimTracks = new TObjArray(); |
222 | ||
223 | ||
224 | fOutputList = new TList(); | |
a62631a9 | 225 | fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) |
30159e6f | 226 | |
f507c09b | 227 | fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data()); |
4d4ce2d2 | 228 | fOADBContainer = new AliOADBContainer("AliEMCALgeo"); |
229 | fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo"); | |
230 | ||
30159e6f | 231 | fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2); |
232 | fOutputList->Add(fEvtSel); | |
0b21f686 | 233 | |
234 | fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5); | |
235 | fOutputList->Add(fNClusEt10); | |
30159e6f | 236 | |
cc57d293 | 237 | fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5); |
238 | fOutputList->Add(fRecoPV); | |
239 | ||
30159e6f | 240 | fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100); |
241 | fOutputList->Add(fPVtxZ); | |
c7bb0b43 | 242 | |
f507c09b | 243 | fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5); |
244 | fOutputList->Add(fTrMultDist); | |
245 | ||
12fbd6e8 | 246 | 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 | 247 | fMCDirPhotonPtEtaPhi->Sumw2(); |
248 | fOutputList->Add(fMCDirPhotonPtEtaPhi); | |
c7bb0b43 | 249 | |
12fbd6e8 | 250 | 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 | 251 | fMCIsoDirPhotonPtEtaPhi->Sumw2(); |
252 | fOutputList->Add(fMCIsoDirPhotonPtEtaPhi); | |
253 | ||
c7bb0b43 | 254 | fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100); |
255 | fDecayPhotonPtMC->Sumw2(); | |
256 | fOutputList->Add(fDecayPhotonPtMC); | |
257 | ||
30159e6f | 258 | 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 | 259 | fOutputList->Add(fCellAbsIdVsAmpl); |
30159e6f | 260 | |
261 | 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); | |
262 | fOutputList->Add(fNClusHighClusE); | |
263 | ||
f9e2362a | 264 | fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",1000,0,100,400,0,4); |
265 | fOutputList->Add(fHigherPtConeM02); | |
266 | ||
f507c09b | 267 | fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100); |
268 | fOutputList->Add(fClusEtMcPt); | |
269 | ||
270 | fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7); | |
271 | fOutputList->Add(fClusMcDetaDphi); | |
272 | ||
273 | fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5); | |
274 | fOutputList->Add(fNClusPerPho); | |
275 | ||
d4f449df | 276 | 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 | 277 | fOutputList->Add(fMcPtInConeBG); |
278 | ||
d4f449df | 279 | 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 | 280 | fOutputList->Add(fMcPtInConeSBG); |
281 | ||
d4f449df | 282 | 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 | 283 | fOutputList->Add(fMcPtInConeBGnoUE); |
284 | ||
d4f449df | 285 | 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 | 286 | fOutputList->Add(fMcPtInConeSBGnoUE); |
287 | ||
81f2660f | 288 | 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);",1000,0,100,600,-10,50); |
289 | fOutputList->Add(fAllIsoEtMcGamma); | |
290 | ||
291 | 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);",1000,0,100,600,-10,50); | |
292 | fOutputList->Add(fAllIsoNoUeEtMcGamma); | |
293 | ||
f912f9a9 | 294 | |
092ceec8 | 295 | fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20); |
296 | fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus); | |
f507c09b | 297 | |
84898d7b | 298 | Int_t nEt=1000, 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 | 299 | Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt}; |
c1f18270 | 300 | fNDimensions = sizeof(bins)/sizeof(Int_t); |
301 | const Int_t ndims = fNDimensions; | |
84898d7b | 302 | Double_t xmin[] = { -0.5, 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-0.5}; |
303 | Double_t xmax[] = { 99.5, 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5}; | |
80c98845 | 304 | if(fPeriod.Contains("11h")){ |
305 | xmax[12]=3999.5; | |
306 | } | |
f507c09b | 307 | 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 | 308 | fOutputList->Add(fHnOutput); |
309 | ||
2b7205ad | 310 | //QA outputs |
311 | fQAList = new TList(); | |
312 | fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging) | |
965c985f | 313 | |
2b7205ad | 314 | fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5); |
315 | fNTracks->Sumw2(); | |
316 | fQAList->Add(fNTracks); | |
317 | ||
318 | fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5); | |
319 | fEmcNCells->Sumw2(); | |
320 | fQAList->Add(fEmcNCells); | |
321 | fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5); | |
322 | fEmcNClus->Sumw2(); | |
323 | fQAList->Add(fEmcNClus); | |
324 | fEmcNClusCut = new TH1F("fEmcNClusCut",";n/event;count",120,-0.5,119.5); | |
325 | fEmcNClusCut->Sumw2(); | |
326 | fQAList->Add(fEmcNClusCut); | |
327 | fNTracksECut = new TH1F("fNTracksECut",";n/event;count",120,-0.5,119.5); | |
328 | fNTracksECut->Sumw2(); | |
329 | fQAList->Add(fNTracksECut); | |
330 | fEmcNCellsCut = new TH1F("fEmcNCellsCut",";n/event;count",120,-0.5,119.5); | |
331 | fEmcNCellsCut->Sumw2(); | |
332 | fQAList->Add(fEmcNCellsCut); | |
e4ba80ad | 333 | fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",100,-0.25,49.75,63,0,6.3); |
334 | fEmcClusEPhi->Sumw2(); | |
335 | fQAList->Add(fEmcClusEPhi); | |
336 | fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",";GeV;#phi",100,-0.25,49.75,63,0,6.3); | |
337 | fEmcClusEPhiCut->Sumw2(); | |
338 | fQAList->Add(fEmcClusEPhiCut); | |
339 | fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",100,-0.25,49.75,19,-0.9,0.9); | |
340 | fEmcClusEEta->Sumw2(); | |
341 | fQAList->Add(fEmcClusEEta); | |
342 | fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",";GeV;#eta",100,-0.25,49.75,18,-0.9,0.9); | |
343 | fEmcClusEEtaCut->Sumw2(); | |
344 | fQAList->Add(fEmcClusEEtaCut); | |
345 | fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3); | |
346 | fTrackPtPhi->Sumw2(); | |
347 | fQAList->Add(fTrackPtPhi); | |
348 | fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3); | |
349 | fTrackPtPhiCut->Sumw2(); | |
350 | fQAList->Add(fTrackPtPhiCut); | |
351 | fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9); | |
352 | fTrackPtEta->Sumw2(); | |
353 | fQAList->Add(fTrackPtEta); | |
354 | fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9); | |
355 | fTrackPtEtaCut->Sumw2(); | |
356 | fQAList->Add(fTrackPtEtaCut); | |
965c985f | 357 | |
3ae97198 | 358 | |
359 | fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",100,-0.25,49.75,63,0,6.3); | |
360 | fMaxCellEPhi->Sumw2(); | |
361 | fQAList->Add(fMaxCellEPhi); | |
362 | ||
30159e6f | 363 | PostData(1, fOutputList); |
2b7205ad | 364 | PostData(2, fQAList); |
30159e6f | 365 | } |
366 | ||
367 | //________________________________________________________________________ | |
368 | void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) | |
369 | { | |
bd092f0f | 370 | // Main loop, called for each event. |
2e20efe5 | 371 | fESDClusters = 0; |
372 | fESDCells = 0; | |
373 | fAODClusters = 0; | |
374 | fAODCells = 0; | |
bd092f0f | 375 | // event trigger selection |
30159e6f | 376 | Bool_t isSelected = 0; |
751194e8 | 377 | if(fPeriod.Contains("11a")){ |
378 | if(fTrigBit.Contains("kEMC")) | |
379 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1); | |
380 | else | |
381 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); | |
382 | } | |
383 | else{ | |
384 | if(fTrigBit.Contains("kEMC")) | |
385 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7); | |
386 | else | |
f507c09b | 387 | if(fTrigBit.Contains("kMB")) |
388 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); | |
389 | else | |
390 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); | |
751194e8 | 391 | } |
80c98845 | 392 | if(fPeriod.Contains("11h")){ |
393 | if(fTrigBit.Contains("kAny")) | |
394 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny); | |
395 | if(fTrigBit.Contains("kAnyINT")) | |
396 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT); | |
397 | } | |
c7bb0b43 | 398 | if(fIsMc) |
399 | isSelected = kTRUE; | |
ecd47673 | 400 | if(fDebug) |
401 | printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc); | |
30159e6f | 402 | if(!isSelected ) |
403 | return; | |
f3843637 | 404 | if(fIsMc){ |
405 | TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree(); | |
406 | TFile *file = (TFile*)tree->GetCurrentFile(); | |
407 | TString filename = file->GetName(); | |
408 | if(!filename.Contains(fPathStrOpt.Data())) | |
409 | return; | |
410 | } | |
2e20efe5 | 411 | AliVEvent *event = (AliVEvent*)InputEvent(); |
412 | if (!event) { | |
413 | printf("ERROR: event not available\n"); | |
30159e6f | 414 | return; |
415 | } | |
822c9944 | 416 | Int_t runnumber = InputEvent()->GetRunNumber() ; |
7a7744db | 417 | if(fDebug) |
418 | printf("run number = %d\n",runnumber); | |
419 | ||
2e20efe5 | 420 | fESD = dynamic_cast<AliESDEvent*>(event); |
7a7744db | 421 | if(!fESD){ |
422 | fAOD = dynamic_cast<AliAODEvent*>(event); | |
423 | if(!fAOD){ | |
424 | printf("ERROR: Invalid type of event!!!\n"); | |
425 | return; | |
426 | } | |
427 | else if(fDebug) | |
428 | printf("AOD EVENT!\n"); | |
429 | } | |
30159e6f | 430 | |
431 | fEvtSel->Fill(0); | |
ecd47673 | 432 | if(fDebug) |
26a8db0f | 433 | printf("event is ok,\n run number=%d\n",runnumber); |
434 | ||
30159e6f | 435 | |
2e20efe5 | 436 | AliVVertex *pv = (AliVVertex*)event->GetPrimaryVertex(); |
437 | Bool_t pvStatus = kTRUE; | |
438 | if(fESD){ | |
439 | AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex(); | |
440 | pvStatus = esdv->GetStatus(); | |
441 | } | |
d1582697 | 442 | if(!pv) |
30159e6f | 443 | return; |
2e20efe5 | 444 | if(!pvStatus) |
d1582697 | 445 | fRecoPV->Fill(0); |
cc57d293 | 446 | else |
447 | fRecoPV->Fill(1); | |
30159e6f | 448 | fPVtxZ->Fill(pv->GetZ()); |
449 | if(TMath::Abs(pv->GetZ())>15) | |
450 | return; | |
ecd47673 | 451 | if(fDebug) |
452 | printf("passed vertex cut\n"); | |
30159e6f | 453 | |
454 | fEvtSel->Fill(1); | |
2e20efe5 | 455 | if(fESD) |
456 | fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks")); | |
457 | if(fAOD) | |
458 | fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks()); | |
bd092f0f | 459 | |
06a28959 | 460 | if(!fTracks){ |
461 | AliError("Track array in event is NULL!"); | |
2e20efe5 | 462 | if(fDebug) |
463 | printf("returning due to not finding tracks!\n"); | |
06a28959 | 464 | return; |
465 | } | |
30159e6f | 466 | // Track loop to fill a pT spectrum |
3f4073ba | 467 | const Int_t Ntracks = fTracks->GetEntriesFast(); |
468 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { | |
469 | // for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { | |
470 | //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks); | |
2e20efe5 | 471 | AliVTrack *track = (AliVTrack*)fTracks->At(iTracks); |
30159e6f | 472 | if (!track) |
473 | continue; | |
2e20efe5 | 474 | AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track); |
475 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
476 | if (esdTrack && fPrTrCuts && fPrTrCuts->IsSelected(track)){ | |
30159e6f | 477 | fSelPrimTracks->Add(track); |
2b7205ad | 478 | /*if(fTrackMaxPt<track->Pt()) |
479 | fTrackMaxPt = track->Pt();*/ | |
30159e6f | 480 | //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi()); |
481 | } | |
2b7205ad | 482 | else if(aodTrack){ |
dcd3dca1 | 483 | if (!aodTrack->IsHybridGlobalConstrainedGlobal()) |
484 | continue ; | |
d8c08ce4 | 485 | /*if(!aodTrack->TestFilterBit(fFilterBit)) |
486 | continue;*/ | |
2e20efe5 | 487 | fSelPrimTracks->Add(track); |
2b7205ad | 488 | /*if(fTrackMaxPt<track->Pt()) |
489 | fTrackMaxPt = track->Pt();*/ | |
490 | } | |
bd092f0f | 491 | } |
492 | ||
4d4ce2d2 | 493 | TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices"); |
822c9944 | 494 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ |
495 | if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) | |
496 | break; | |
497 | /*if(fESD) | |
498 | fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod); | |
499 | else*/ | |
500 | // if(event->GetEMCALMatrix(mod)) | |
501 | fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod); | |
502 | fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod); | |
30159e6f | 503 | } |
4d4ce2d2 | 504 | |
2e20efe5 | 505 | if(fESD){ |
506 | AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); | |
507 | fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); | |
508 | if(fDebug) | |
509 | printf("ESD Track mult= %d\n",fTrackMult); | |
510 | } | |
511 | else if(fAOD){ | |
512 | fTrackMult = Ntracks; | |
513 | if(fDebug) | |
514 | printf("AOD Track mult= %d\n",fTrackMult); | |
515 | } | |
f507c09b | 516 | fTrMultDist->Fill(fTrackMult); |
16a4050e | 517 | |
2e20efe5 | 518 | if(fESD){ |
519 | TList *l = fESD->GetList(); | |
520 | fESDClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
521 | fESDCells = fESD->GetEMCALCells(); | |
522 | if(fDebug) | |
523 | printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast()); | |
524 | } | |
525 | else if(fAOD){ | |
526 | fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters()); | |
527 | fAODCells = fAOD->GetEMCALCells(); | |
528 | if(fDebug) | |
529 | printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast()); | |
530 | } | |
30159e6f | 531 | |
c7bb0b43 | 532 | |
533 | fMCEvent = MCEvent(); | |
ecd47673 | 534 | if(fMCEvent){ |
535 | if(fDebug) | |
ec9ab86b | 536 | std::cout<<"MCevent exists\n"; |
c7bb0b43 | 537 | fStack = (AliStack*)fMCEvent->Stack(); |
ecd47673 | 538 | } |
539 | else{ | |
540 | if(fDebug) | |
ec9ab86b | 541 | std::cout<<"ERROR: NO MC EVENT!!!!!!\n"; |
ecd47673 | 542 | } |
112eb594 | 543 | LoopOnCells(); |
f507c09b | 544 | FollowGamma(); |
545 | if(fDebug) | |
546 | printf("passed calling of FollowGamma\n"); | |
547 | FillClusHists(); | |
548 | if(fDebug) | |
549 | printf("passed calling of FillClusHists\n"); | |
63e40cb8 | 550 | FillMcHists(); |
ecd47673 | 551 | if(fDebug) |
552 | printf("passed calling of FillMcHists\n"); | |
ab031886 | 553 | if(fESD) |
554 | FillQA(); | |
2b7205ad | 555 | if(fDebug) |
556 | printf("passed calling of FillQA\n"); | |
2e20efe5 | 557 | /*if(fESD) |
558 | fESDClusters->Clear();*/ | |
f507c09b | 559 | fSelPrimTracks->Clear(); |
560 | fNClusForDirPho = 0; | |
112eb594 | 561 | fNCells50 = 0; |
c7bb0b43 | 562 | |
30159e6f | 563 | PostData(1, fOutputList); |
2b7205ad | 564 | PostData(2, fQAList); |
30159e6f | 565 | } |
bd092f0f | 566 | |
30159e6f | 567 | //________________________________________________________________________ |
568 | void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() | |
569 | { | |
2e20efe5 | 570 | if(fDebug) |
571 | printf("Inside FillClusHists()....\n"); | |
bd092f0f | 572 | // Fill cluster histograms. |
2e20efe5 | 573 | TObjArray *clusters = fESDClusters; |
bd092f0f | 574 | |
2e20efe5 | 575 | if (!clusters){ |
576 | clusters = fAODClusters; | |
577 | if(fDebug) | |
578 | printf("ESD clusters empty..."); | |
579 | } | |
580 | if (!clusters){ | |
581 | if(fDebug) | |
582 | printf(" and AOD clusters as well!!!\n"); | |
30159e6f | 583 | return; |
2e20efe5 | 584 | } |
585 | if(fDebug) | |
586 | printf("\n"); | |
587 | ||
588 | const Int_t nclus = clusters->GetEntries(); | |
30159e6f | 589 | if(nclus==0) |
590 | return; | |
ecd47673 | 591 | if(fDebug) |
592 | printf("Inside FillClusHists and there are %d clusters\n",nclus); | |
e681d9ce | 593 | Double_t maxE = 0; |
0b21f686 | 594 | Int_t nclus10 = 0; |
f507c09b | 595 | Double_t ptmc=-1; |
30159e6f | 596 | for(Int_t ic=0;ic<nclus;ic++){ |
597 | maxE=0; | |
2e20efe5 | 598 | AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic)); |
30159e6f | 599 | if(!c) |
600 | continue; | |
601 | if(!c->IsEMCAL()) | |
602 | continue; | |
965c985f | 603 | if(c->E()<fECut) |
604 | continue; | |
30159e6f | 605 | Short_t id; |
606 | Double_t Emax = GetMaxCellEnergy( c, id); | |
607 | Double_t Ecross = GetCrossEnergy( c, id); | |
608 | if((1-Ecross/Emax)>fExoticCut) | |
609 | continue; | |
610 | Float_t clsPos[3] = {0,0,0}; | |
611 | c->GetPosition(clsPos); | |
612 | TVector3 clsVec(clsPos); | |
3dc2825e | 613 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); |
f507c09b | 614 | if(fDebug) |
615 | printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E()); | |
0b21f686 | 616 | if(Et>10) |
617 | nclus10++; | |
30159e6f | 618 | Float_t ceiso, cephiband, cecore; |
619 | Float_t triso, trphiband, trcore; | |
620 | Float_t alliso, allphiband, allcore; | |
621 | Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR; | |
622 | Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04); | |
f224025f | 623 | GetCeIso(clsVec, id, ceiso, cephiband, cecore); |
30159e6f | 624 | GetTrIso(clsVec, triso, trphiband, trcore); |
f9e2362a | 625 | Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz()); |
626 | if(Et>10 && Et<15 && dr>0.025){ | |
627 | fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02()); | |
628 | if(fDebug) | |
629 | printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02()); | |
630 | } | |
30159e6f | 631 | alliso = ceiso + triso; |
632 | allphiband = cephiband + trphiband; | |
633 | allcore = cecore + trcore; | |
634 | Float_t ceisoue = cephiband/phibandArea*netConeArea; | |
635 | Float_t trisoue = trphiband/phibandArea*netConeArea; | |
636 | Float_t allisoue = allphiband/phibandArea*netConeArea; | |
f912f9a9 | 637 | Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR); |
638 | if(fDebug && Et>10) | |
639 | printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et); | |
640 | if(c->GetM02()>0.5 && c->GetM02()<2.0){ | |
641 | fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum); | |
642 | fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum); | |
643 | } | |
b16ded35 | 644 | if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){ |
f912f9a9 | 645 | fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum); |
646 | fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum); | |
81f2660f | 647 | if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){ |
672df96f | 648 | fAllIsoEtMcGamma->Fill(Et, alliso-/*Et*/cecore-allisoue); |
649 | fAllIsoNoUeEtMcGamma->Fill(Et, alliso-/*Et*/cecore); | |
81f2660f | 650 | } |
f912f9a9 | 651 | } |
965c985f | 652 | const Int_t ndims = fNDimensions; |
653 | Double_t outputValues[ndims]; | |
e53ab710 | 654 | if(mcptsum<2) |
655 | ptmc = GetClusSource(c); | |
656 | else | |
657 | ptmc = -0.1; | |
965c985f | 658 | outputValues[0] = Et; |
659 | outputValues[1] = c->GetM02(); | |
1c86c72c | 660 | outputValues[2] = ceiso/*cecore*/-ceisoue; |
965c985f | 661 | outputValues[3] = triso-trisoue; |
1c86c72c | 662 | outputValues[4] = alliso/*cecore*/-allisoue - trcore; |
663 | outputValues[5] = ceiso; | |
664 | outputValues[6] = alliso - trcore; | |
717bb45b | 665 | outputValues[7] = c->GetTrackDx(); |
666 | outputValues[8] = c->GetTrackDz(); | |
667 | outputValues[9] = clsVec.Eta(); | |
668 | outputValues[10] = clsVec.Phi(); | |
2e20efe5 | 669 | if(fESDCells) |
670 | outputValues[11] = fESDCells->GetCellTime(id); | |
671 | else if(fAODCells) | |
672 | outputValues[11] = fAODCells->GetCellTime(id); | |
16a4050e | 673 | outputValues[12] = fTrackMult; |
f507c09b | 674 | outputValues[13] = ptmc; |
965c985f | 675 | fHnOutput->Fill(outputValues); |
30159e6f | 676 | if(c->E()>maxE) |
677 | maxE = c->E(); | |
678 | } | |
679 | fNClusHighClusE->Fill(maxE,nclus); | |
ee3d9c10 | 680 | fMaxEClus = maxE; |
0b21f686 | 681 | fNClusEt10->Fill(nclus10); |
f507c09b | 682 | fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho); |
30159e6f | 683 | } |
bd092f0f | 684 | |
30159e6f | 685 | //________________________________________________________________________ |
f224025f | 686 | void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core) |
30159e6f | 687 | { |
bd092f0f | 688 | // Get cell isolation. |
2e20efe5 | 689 | AliVCaloCells *cells = fESDCells; |
1c86c72c | 690 | if (!cells){ |
2e20efe5 | 691 | cells = fAODCells; |
1c86c72c | 692 | if(fDebug) |
693 | printf("ESD cells empty..."); | |
694 | } | |
695 | if (!cells){ | |
696 | if(fDebug) | |
697 | printf(" and AOD cells are empty as well!!!\n"); | |
30159e6f | 698 | return; |
1c86c72c | 699 | } |
2e20efe5 | 700 | |
1c86c72c | 701 | TObjArray *clusters = fESDClusters; |
702 | if (!clusters) | |
703 | clusters = fAODClusters; | |
704 | if (!clusters) | |
705 | return; | |
706 | ||
707 | ||
708 | const Int_t nclus = clusters->GetEntries(); | |
709 | //const Int_t ncells = cells->GetNumberOfCells(); | |
30159e6f | 710 | Float_t totiso=0; |
711 | Float_t totband=0; | |
712 | Float_t totcore=0; | |
713 | Float_t etacl = vec.Eta(); | |
714 | Float_t phicl = vec.Phi(); | |
2e20efe5 | 715 | Float_t maxtcl = cells->GetCellTime(maxid); |
30159e6f | 716 | if(phicl<0) |
717 | phicl+=TMath::TwoPi(); | |
1c86c72c | 718 | /*Int_t absid = -1; |
30159e6f | 719 | Float_t eta=-1, phi=-1; |
720 | for(int icell=0;icell<ncells;icell++){ | |
2e20efe5 | 721 | absid = TMath::Abs(cells->GetCellNumber(icell)); |
722 | Float_t celltime = cells->GetCellTime(absid); | |
f224025f | 723 | //if(TMath::Abs(celltime)>2e-8 && (!fIsMc)) |
724 | if(TMath::Abs(celltime-maxtcl)>2e-8 ) | |
c3f33463 | 725 | continue; |
30159e6f | 726 | if(!fGeom) |
727 | return; | |
728 | fGeom->EtaPhiFromIndex(absid,eta,phi); | |
729 | Float_t dphi = TMath::Abs(phi-phicl); | |
730 | Float_t deta = TMath::Abs(eta-etacl); | |
1c86c72c | 731 | Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/ |
732 | for(int ic=0;ic<nclus;ic++){ | |
733 | AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic)); | |
734 | if(!c) | |
735 | continue; | |
736 | if(!c->IsEMCAL()) | |
737 | continue; | |
738 | Short_t id; | |
2280356c | 739 | GetMaxCellEnergy( c, id); |
1c86c72c | 740 | Double_t maxct = cells->GetCellTime(id); |
741 | if(TMath::Abs(maxtcl-maxct)>2.5e-9) | |
742 | continue; | |
743 | Float_t clsPos[3] = {0,0,0}; | |
744 | c->GetPosition(clsPos); | |
745 | TVector3 cv(clsPos); | |
746 | Double_t Et = c->E()*TMath::Sin(cv.Theta()); | |
747 | Float_t dphi = TMath::Abs(cv.Phi()-phicl); | |
748 | Float_t deta = TMath::Abs(cv.Eta()-etacl); | |
30159e6f | 749 | Float_t R = TMath::Sqrt(deta*deta + dphi*dphi); |
1c86c72c | 750 | if(R<0.007) |
751 | continue; | |
22814624 | 752 | if(maxid==id) |
753 | continue; | |
1c86c72c | 754 | Double_t matchedpt = GetTrackMatchedPt(c->GetTrackMatchedIndex()); |
755 | Double_t nEt = TMath::Max(Et-matchedpt, 0.0); | |
756 | if(nEt<0) | |
757 | printf("nEt=%1.1f\n",nEt); | |
30159e6f | 758 | if(R<fIsoConeR){ |
1c86c72c | 759 | totiso += nEt; |
30159e6f | 760 | if(R<0.04) |
1c86c72c | 761 | totcore += nEt; |
30159e6f | 762 | } |
763 | else{ | |
764 | if(dphi>fIsoConeR) | |
765 | continue; | |
766 | if(deta<fIsoConeR) | |
767 | continue; | |
1c86c72c | 768 | totband += nEt; |
30159e6f | 769 | } |
770 | } | |
771 | iso = totiso; | |
772 | phiband = totband; | |
773 | core = totcore; | |
774 | } | |
775 | //________________________________________________________________________ | |
776 | void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core) | |
777 | { | |
bd092f0f | 778 | // Get track isolation. |
779 | ||
30159e6f | 780 | if(!fSelPrimTracks) |
781 | return; | |
f9e2362a | 782 | fHigherPtCone = 0; |
30159e6f | 783 | const Int_t ntracks = fSelPrimTracks->GetEntries(); |
784 | Double_t totiso=0; | |
785 | Double_t totband=0; | |
786 | Double_t totcore=0; | |
787 | Double_t etacl = vec.Eta(); | |
788 | Double_t phicl = vec.Phi(); | |
789 | if(phicl<0) | |
790 | phicl+=TMath::TwoPi(); | |
791 | for(int itrack=0;itrack<ntracks;itrack++){ | |
3f4073ba | 792 | AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack)); |
30159e6f | 793 | if(!track) |
794 | continue; | |
795 | Double_t dphi = TMath::Abs(track->Phi()-phicl); | |
796 | Double_t deta = TMath::Abs(track->Eta()-etacl); | |
797 | Double_t R = TMath::Sqrt(deta*deta + dphi*dphi); | |
798 | Double_t pt = track->Pt(); | |
f9e2362a | 799 | if(pt>fHigherPtCone) |
800 | fHigherPtCone = pt; | |
30159e6f | 801 | if(R<fIsoConeR){ |
802 | totiso += track->Pt(); | |
803 | if(R<0.04) | |
804 | totcore += pt; | |
805 | } | |
806 | else{ | |
807 | if(dphi>fIsoConeR) | |
808 | continue; | |
809 | if(deta<fIsoConeR) | |
810 | continue; | |
811 | totband += track->Pt(); | |
812 | } | |
813 | } | |
814 | iso = totiso; | |
815 | phiband = totband; | |
816 | core = totcore; | |
817 | } | |
bd092f0f | 818 | |
30159e6f | 819 | //________________________________________________________________________ |
820 | Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) | |
821 | { | |
822 | // Calculate the energy of cross cells around the leading cell. | |
823 | ||
824 | AliVCaloCells *cells = 0; | |
2e20efe5 | 825 | cells = fESDCells; |
826 | if (!cells) | |
827 | cells = fAODCells; | |
30159e6f | 828 | if (!cells) |
829 | return 0; | |
830 | ||
30159e6f | 831 | if (!fGeom) |
832 | return 0; | |
833 | ||
834 | Int_t iSupMod = -1; | |
835 | Int_t iTower = -1; | |
836 | Int_t iIphi = -1; | |
837 | Int_t iIeta = -1; | |
838 | Int_t iphi = -1; | |
839 | Int_t ieta = -1; | |
840 | Int_t iphis = -1; | |
841 | Int_t ietas = -1; | |
842 | ||
843 | Double_t crossEnergy = 0; | |
844 | ||
845 | fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta); | |
846 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas); | |
847 | ||
848 | Int_t ncells = cluster->GetNCells(); | |
849 | for (Int_t i=0; i<ncells; i++) { | |
850 | Int_t cellAbsId = cluster->GetCellAbsId(i); | |
851 | fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); | |
852 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); | |
853 | Int_t aphidiff = TMath::Abs(iphi-iphis); | |
854 | if (aphidiff>1) | |
855 | continue; | |
856 | Int_t aetadiff = TMath::Abs(ieta-ietas); | |
857 | if (aetadiff>1) | |
858 | continue; | |
859 | if ( (aphidiff==1 && aetadiff==0) || | |
860 | (aphidiff==0 && aetadiff==1) ) { | |
861 | crossEnergy += cells->GetCellAmplitude(cellAbsId); | |
862 | } | |
863 | } | |
864 | ||
865 | return crossEnergy; | |
866 | } | |
867 | ||
30159e6f | 868 | //________________________________________________________________________ |
869 | Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const | |
870 | { | |
871 | // Get maximum energy of attached cell. | |
872 | ||
873 | id = -1; | |
874 | ||
875 | AliVCaloCells *cells = 0; | |
2e20efe5 | 876 | cells = fESDCells; |
30159e6f | 877 | if (!cells) |
2e20efe5 | 878 | cells = fAODCells; |
879 | if(!cells) | |
30159e6f | 880 | return 0; |
881 | ||
882 | Double_t maxe = 0; | |
883 | Int_t ncells = cluster->GetNCells(); | |
884 | for (Int_t i=0; i<ncells; i++) { | |
885 | Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); | |
886 | if (e>maxe) { | |
887 | maxe = e; | |
888 | id = cluster->GetCellAbsId(i); | |
889 | } | |
890 | } | |
891 | return maxe; | |
892 | } | |
bd092f0f | 893 | |
c7bb0b43 | 894 | //________________________________________________________________________ |
895 | void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists() | |
896 | { | |
897 | if(!fStack) | |
898 | return; | |
899 | Int_t nTracks = fStack->GetNtrack(); | |
ecd47673 | 900 | if(fDebug) |
901 | printf("Inside FillMcHists and there are %d mcparts\n",nTracks); | |
c7bb0b43 | 902 | for(Int_t iTrack=0;iTrack<nTracks;iTrack++){ |
903 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); | |
904 | if(!mcp) | |
905 | continue; | |
906 | Int_t pdg = mcp->GetPdgCode(); | |
907 | if(pdg!=22) | |
908 | continue; | |
f2acdbe9 | 909 | if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2) |
910 | continue; | |
c7bb0b43 | 911 | Int_t imom = mcp->GetMother(0); |
912 | if(imom<0 || imom>nTracks) | |
913 | continue; | |
914 | TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom)); | |
915 | if(!mcmom) | |
916 | continue; | |
917 | Int_t pdgMom = mcmom->GetPdgCode(); | |
ec9ab86b | 918 | if((imom==6 || imom==7) && pdgMom==22) { |
caaf99d3 | 919 | fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
ca03dc51 | 920 | Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR); |
e53ab710 | 921 | if(ptsum<2) |
922 | fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); | |
f507c09b | 923 | if(fNClusForDirPho==0) |
092ceec8 | 924 | fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
40c4f1bf | 925 | if(fDebug){ |
f2acdbe9 | 926 | 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()); |
40c4f1bf | 927 | 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()); |
928 | } | |
f2acdbe9 | 929 | } |
c7bb0b43 | 930 | else{ |
931 | if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000) | |
932 | fDecayPhotonPtMC->Fill(mcp->Pt()); | |
933 | } | |
934 | } | |
935 | } | |
30159e6f | 936 | //________________________________________________________________________ |
f507c09b | 937 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) |
938 | { | |
939 | if(!c) | |
84898d7b | 940 | return -0.1; |
f507c09b | 941 | if(!fStack) |
84898d7b | 942 | return -0.1; |
f507c09b | 943 | Int_t nmcp = fStack->GetNtrack(); |
944 | Int_t clabel = c->GetLabel(); | |
945 | if(fDebug && fMcIdFamily.Contains(Form("%d",clabel))) | |
946 | printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel); | |
947 | if(!fMcIdFamily.Contains(Form("%d",clabel))) | |
84898d7b | 948 | return -0.1; |
f507c09b | 949 | fNClusForDirPho++; |
950 | TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1); | |
951 | Int_t partonpos = partonposstr.Atoi(); | |
952 | if(fDebug) | |
953 | printf("\t^^^^ parton position = %d ^^^^\n",partonpos); | |
954 | if(clabel<0 || clabel>nmcp) | |
2cd24808 | 955 | return -0.1; |
f507c09b | 956 | Float_t clsPos[3] = {0,0,0}; |
957 | c->GetPosition(clsPos); | |
958 | TVector3 clsVec(clsPos); | |
959 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); | |
960 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos)); | |
961 | if(!mcp) | |
84898d7b | 962 | return -0.1; |
f507c09b | 963 | if(fDebug){ |
964 | 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); | |
965 | } | |
966 | Int_t lab1 = mcp->GetFirstDaughter(); | |
967 | if(lab1<0 || lab1>nmcp) | |
84898d7b | 968 | return -0.1; |
f507c09b | 969 | TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1)); |
970 | if(!mcd) | |
84898d7b | 971 | return -0.1; |
f507c09b | 972 | if(fDebug) |
973 | 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); | |
974 | if(mcd->GetPdgCode()==22){ | |
975 | fClusEtMcPt->Fill(mcd->Pt(), Et); | |
976 | fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi()); | |
977 | } | |
978 | else{ | |
979 | printf("Warning: daughter of photon parton is not a photon\n"); | |
84898d7b | 980 | return -0.1; |
f507c09b | 981 | } |
982 | return fDirPhoPt; | |
983 | } | |
984 | //________________________________________________________________________ | |
985 | void AliAnalysisTaskEMCALIsoPhoton::FollowGamma() | |
986 | { | |
987 | if(!fStack) | |
988 | return; | |
989 | Int_t selfid = 6; | |
990 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid)); | |
991 | if(!mcp) | |
992 | return; | |
993 | if(mcp->GetPdgCode()!=22){ | |
994 | mcp = static_cast<TParticle*>(fStack->Particle(++selfid)); | |
995 | if(!mcp) | |
996 | return; | |
997 | } | |
998 | Int_t daug0f = mcp->GetFirstDaughter(); | |
999 | Int_t daug0l = mcp->GetLastDaughter(); | |
1000 | Int_t nd0 = daug0l - daug0f; | |
1001 | if(fDebug) | |
1002 | 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); | |
1003 | fMcIdFamily = Form("%d,",selfid); | |
1004 | GetDaughtersInfo(daug0f,daug0l, selfid,""); | |
1005 | if(fDebug){ | |
1006 | printf("\t---- end of the prompt gamma's family tree ----\n\n"); | |
1007 | printf("Family id string = %s,\n\n",fMcIdFamily.Data()); | |
1008 | } | |
1009 | TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f)); | |
1010 | if(!md) | |
1011 | return; | |
1012 | fDirPhoPt = md->Pt(); | |
1013 | } | |
1014 | //________________________________________________________________________ | |
22ad7981 | 1015 | void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind) |
f507c09b | 1016 | { |
1017 | int nmcp = fStack->GetNtrack(); | |
1018 | if(firstd<0 || firstd>nmcp) | |
1019 | return; | |
1020 | if(lastd<0 || lastd>nmcp) | |
1021 | return; | |
1022 | if(firstd>lastd){ | |
1023 | printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd); | |
1024 | return; | |
1025 | } | |
1026 | TString indenter = Form("\t%s",inputind); | |
1027 | TParticle *dp = 0x0; | |
1028 | if(fDebug) | |
1029 | printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid); | |
1030 | for(int id=firstd; id<lastd+1; id++){ | |
1031 | dp = static_cast<TParticle*>(fStack->Particle(id)); | |
1032 | if(!dp) | |
1033 | continue; | |
1034 | Int_t dfd = dp->GetFirstDaughter(); | |
1035 | Int_t dld = dp->GetLastDaughter(); | |
1036 | Int_t dnd = dld - dfd + 1; | |
1037 | Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy()); | |
1038 | if(fDebug) | |
1039 | 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); | |
1040 | fMcIdFamily += Form("%d,",id); | |
1041 | GetDaughtersInfo(dfd,dld,id,indenter.Data()); | |
1042 | } | |
1043 | } | |
f912f9a9 | 1044 | |
1045 | //________________________________________________________________________ | |
1046 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R) | |
1047 | { | |
1048 | if(!fStack) | |
1049 | return 0; | |
1050 | if(fDebug) | |
1051 | printf("Inside GetMcPtSumInCone!!\n"); | |
1052 | Int_t nTracks = fStack->GetNtrack(); | |
1053 | Float_t ptsum = 0; | |
672df96f | 1054 | TString addpartlabels = ""; |
e53ab710 | 1055 | for(Int_t iTrack=8;iTrack<nTracks;iTrack++){ |
f912f9a9 | 1056 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); |
1057 | if(!mcp) | |
1058 | continue; | |
e53ab710 | 1059 | Int_t status = mcp->GetStatusCode(); |
1060 | if(status!=1) | |
1061 | continue; | |
1062 | Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz()); | |
1063 | if(mcvr>1) | |
f912f9a9 | 1064 | continue; |
1065 | else { | |
1066 | if(fDebug) | |
1067 | printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz()); | |
1068 | } | |
1069 | Float_t dphi = mcp->Phi() - phiclus; | |
1070 | Float_t deta = mcp->Eta() - etaclus; | |
1071 | if(fDebug) | |
1072 | printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta); | |
e53ab710 | 1073 | if(deta>R || dphi>R) |
f912f9a9 | 1074 | continue; |
1075 | Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta); | |
1076 | if(dR>R) | |
1077 | continue; | |
672df96f | 1078 | addpartlabels += Form("%d,",iTrack); |
e53ab710 | 1079 | if(mcp->Pt()<0.2) |
1080 | continue; | |
f912f9a9 | 1081 | ptsum += mcp->Pt(); |
1082 | } | |
1083 | return ptsum; | |
1084 | } | |
f507c09b | 1085 | //________________________________________________________________________ |
2b7205ad | 1086 | void AliAnalysisTaskEMCALIsoPhoton::FillQA() |
1087 | { | |
1088 | if(!fSelPrimTracks) | |
1089 | return; | |
1090 | const int ntracks = fSelPrimTracks->GetEntriesFast(); | |
112eb594 | 1091 | const int ncells = fNCells50;//fESDCells->GetNumberOfCells(); |
2b7205ad | 1092 | const Int_t nclus = fESDClusters->GetEntries(); |
1093 | ||
1094 | fNTracks->Fill(ntracks); | |
1095 | fEmcNCells->Fill(ncells); | |
1096 | fEmcNClus->Fill(nclus); | |
1097 | if(fMaxEClus>fECut){ | |
1098 | fNTracksECut->Fill(ntracks); | |
1099 | fEmcNCellsCut->Fill(ncells); | |
1100 | fEmcNClusCut->Fill(nclus); | |
1101 | } | |
1102 | for(int it=0;it<ntracks;it++){ | |
1103 | AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it); | |
1104 | if(!t) | |
1105 | continue; | |
e4ba80ad | 1106 | fTrackPtPhi->Fill(t->Pt(),t->Phi()); |
1107 | fTrackPtEta->Fill(t->Pt(),t->Eta()); | |
1108 | if(fMaxEClus>fECut){ | |
1109 | fTrackPtPhiCut->Fill(t->Pt(), t->Phi()); | |
1110 | fTrackPtEtaCut->Fill(t->Pt(), t->Eta()); | |
1111 | } | |
2b7205ad | 1112 | } |
1113 | for(int ic=0;ic<nclus;ic++){ | |
3ae97198 | 1114 | AliVCluster *c = dynamic_cast<AliVCluster*>(fESDClusters->At(ic)); |
1115 | //AliESDCaloCluster *c = (AliESDCaloCluster*)fESDClusters->At(ic); | |
2b7205ad | 1116 | if(!c) |
1117 | continue; | |
3ae97198 | 1118 | if(!c->IsEMCAL()) |
1119 | continue; | |
e4ba80ad | 1120 | Float_t clsPos[3] = {0,0,0}; |
1121 | c->GetPosition(clsPos); | |
1122 | TVector3 clsVec(clsPos); | |
1123 | Double_t cphi = clsVec.Phi(); | |
1124 | Double_t ceta = clsVec.Eta(); | |
1125 | fEmcClusEPhi->Fill(c->E(), cphi); | |
1126 | fEmcClusEEta->Fill(c->E(), ceta); | |
1127 | if(fMaxEClus>fECut){ | |
1128 | fEmcClusEPhiCut->Fill(c->E(), cphi); | |
1129 | fEmcClusEEtaCut->Fill(c->E(), ceta); | |
1130 | } | |
2b7205ad | 1131 | } |
1132 | } | |
1133 | //________________________________________________________________________ | |
1c86c72c | 1134 | Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex) |
1135 | { | |
1136 | Double_t pt = 0; | |
1137 | if(!fTracks) | |
1138 | return pt; | |
1139 | if(matchIndex<0 || matchIndex>fTracks->GetEntries()){ | |
1140 | if(fDebug) | |
1141 | printf("track-matched index out of track array range!!!\n"); | |
1142 | return pt; | |
1143 | } | |
1144 | AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex)); | |
1145 | if(!track){ | |
1146 | if(fDebug) | |
1147 | printf("track-matched pointer does not exist!!!\n"); | |
1148 | return pt; | |
1149 | } | |
1150 | if(fESD){ | |
1151 | if(!fPrTrCuts) | |
1152 | return pt; | |
1153 | if(!fPrTrCuts->IsSelected(track)) | |
1154 | return pt; | |
1155 | pt = track->Pt(); | |
1156 | } | |
1157 | return pt; | |
1158 | } | |
1159 | //________________________________________________________________________ | |
112eb594 | 1160 | void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells() |
1161 | { | |
1162 | AliVCaloCells *cells = 0; | |
1163 | cells = fESDCells; | |
1164 | if (!cells) | |
1165 | cells = fAODCells; | |
1166 | if(!cells) | |
1167 | return; | |
3ae97198 | 1168 | Double_t maxe = 0; |
1169 | Double_t maxphi = -10; | |
112eb594 | 1170 | Int_t ncells = cells->GetNumberOfCells(); |
3ae97198 | 1171 | Double_t eta,phi; |
112eb594 | 1172 | for (Int_t i=0; i<ncells; i++) { |
1173 | Short_t absid = TMath::Abs(cells->GetCellNumber(i)); | |
1174 | Double_t e = cells->GetCellAmplitude(absid); | |
1175 | if(e>0.05) | |
1176 | fNCells50++; | |
3ae97198 | 1177 | else |
1178 | continue; | |
1179 | fGeom->EtaPhiFromIndex(absid,eta,phi); | |
1180 | if(maxe<e){ | |
1181 | maxe = e; | |
1182 | maxphi = phi; | |
1183 | } | |
112eb594 | 1184 | } |
3ae97198 | 1185 | fMaxCellEPhi->Fill(maxe,maxphi); |
112eb594 | 1186 | |
1187 | } | |
1188 | //________________________________________________________________________ | |
30159e6f | 1189 | void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) |
1190 | { | |
bd092f0f | 1191 | // Called once at the end of the query. |
30159e6f | 1192 | } |