]>
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" | |
bd092f0f | 27 | #include "AliMCEvent.h" |
28 | #include "AliMCEventHandler.h" | |
29 | #include "AliStack.h" | |
30159e6f | 30 | #include "AliV0vertexer.h" |
30159e6f | 31 | #include "AliVCluster.h" |
32 | ||
30159e6f | 33 | ClassImp(AliAnalysisTaskEMCALIsoPhoton) |
34 | ||
35 | //________________________________________________________________________ | |
bd092f0f | 36 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : |
37 | AliAnalysisTaskSE(), | |
38 | fCaloClusters(0), | |
39 | fSelPrimTracks(0), | |
3f4073ba | 40 | fTracks(0), |
bd092f0f | 41 | fEMCalCells(0), |
42 | fPrTrCuts(0), | |
43 | fGeom(0x0), | |
9148fa89 | 44 | fGeoName("EMCAL_COMPLETEV1"), |
45 | fPeriod("LHC11c"), | |
46 | fTrigBit("kEMC7"), | |
47 | fIsTrain(0), | |
c7bb0b43 | 48 | fIsMc(0), |
ecd47673 | 49 | fDebug(0), |
f3843637 | 50 | fPathStrOpt("/"), |
9148fa89 | 51 | fExoticCut(0.97), |
52 | fIsoConeR(0.4), | |
965c985f | 53 | fNDimensions(7), |
54 | fECut(3.), | |
16a4050e | 55 | fTrackMult(0), |
f507c09b | 56 | fMcIdFamily(""), |
57 | fNClusForDirPho(0), | |
58 | fDirPhoPt(0), | |
f9e2362a | 59 | fHigherPtCone(0), |
bd092f0f | 60 | fESD(0), |
c7bb0b43 | 61 | fMCEvent(0), |
62 | fStack(0), | |
bd092f0f | 63 | fOutputList(0), |
64 | fEvtSel(0), | |
0b21f686 | 65 | fNClusEt10(0), |
cc57d293 | 66 | fRecoPV(0), |
f507c09b | 67 | fPVtxZ(0), |
68 | fTrMultDist(0), | |
caaf99d3 | 69 | fMCDirPhotonPtEtaPhi(0), |
c7bb0b43 | 70 | fDecayPhotonPtMC(0), |
bd092f0f | 71 | fCellAbsIdVsAmpl(0), |
16a4050e | 72 | fNClusHighClusE(0), |
f9e2362a | 73 | fHigherPtConeM02(0), |
f507c09b | 74 | fClusEtMcPt(0), |
75 | fClusMcDetaDphi(0), | |
76 | fNClusPerPho(0), | |
f912f9a9 | 77 | fMcPtInConeBG(0), |
78 | fMcPtInConeSBG(0), | |
79 | fMcPtInConeBGnoUE(0), | |
80 | fMcPtInConeSBGnoUE(0), | |
092ceec8 | 81 | fMCDirPhotonPtEtaPhiNoClus(0), |
965c985f | 82 | fHnOutput(0) |
bd092f0f | 83 | { |
84 | // Default constructor. | |
85 | } | |
86 | ||
87 | //________________________________________________________________________ | |
88 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : | |
89 | AliAnalysisTaskSE(name), | |
30159e6f | 90 | fCaloClusters(0), |
91 | fSelPrimTracks(0), | |
3f4073ba | 92 | fTracks(0), |
30159e6f | 93 | fEMCalCells(0), |
94 | fPrTrCuts(0), | |
95 | fGeom(0x0), | |
96 | fGeoName("EMCAL_COMPLETEV1"), | |
97 | fPeriod("LHC11c"), | |
751194e8 | 98 | fTrigBit("kEMC7"), |
30159e6f | 99 | fIsTrain(0), |
c7bb0b43 | 100 | fIsMc(0), |
ecd47673 | 101 | fDebug(0), |
f3843637 | 102 | fPathStrOpt("/"), |
30159e6f | 103 | fExoticCut(0.97), |
104 | fIsoConeR(0.4), | |
965c985f | 105 | fNDimensions(7), |
106 | fECut(3.), | |
16a4050e | 107 | fTrackMult(0), |
f507c09b | 108 | fMcIdFamily(""), |
109 | fNClusForDirPho(0), | |
110 | fDirPhoPt(0), | |
f9e2362a | 111 | fHigherPtCone(0), |
30159e6f | 112 | fESD(0), |
c7bb0b43 | 113 | fMCEvent(0), |
114 | fStack(0), | |
30159e6f | 115 | fOutputList(0), |
30159e6f | 116 | fEvtSel(0), |
0b21f686 | 117 | fNClusEt10(0), |
cc57d293 | 118 | fRecoPV(0), |
bd092f0f | 119 | fPVtxZ(0), |
f507c09b | 120 | fTrMultDist(0), |
caaf99d3 | 121 | fMCDirPhotonPtEtaPhi(0), |
c7bb0b43 | 122 | fDecayPhotonPtMC(0), |
123 | fCellAbsIdVsAmpl(0), | |
bd092f0f | 124 | fNClusHighClusE(0), |
f9e2362a | 125 | fHigherPtConeM02(0), |
f507c09b | 126 | fClusEtMcPt(0), |
127 | fClusMcDetaDphi(0), | |
128 | fNClusPerPho(0), | |
f912f9a9 | 129 | fMcPtInConeBG(0), |
130 | fMcPtInConeSBG(0), | |
131 | fMcPtInConeBGnoUE(0), | |
132 | fMcPtInConeSBGnoUE(0), | |
092ceec8 | 133 | fMCDirPhotonPtEtaPhiNoClus(0), |
965c985f | 134 | fHnOutput(0) |
30159e6f | 135 | { |
136 | // Constructor | |
137 | ||
138 | // Define input and output slots here | |
139 | // Input slot #0 works with a TChain | |
140 | DefineInput(0, TChain::Class()); | |
141 | // Output slot #0 id reserved by the base class for AOD | |
142 | // Output slot #1 writes into a TH1 container | |
143 | DefineOutput(1, TList::Class()); | |
144 | } | |
bd092f0f | 145 | |
30159e6f | 146 | //________________________________________________________________________ |
147 | void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() | |
148 | { | |
bd092f0f | 149 | // Create histograms, called once. |
30159e6f | 150 | |
151 | fCaloClusters = new TRefArray(); | |
152 | fSelPrimTracks = new TObjArray(); | |
153 | ||
154 | ||
155 | fOutputList = new TList(); | |
a62631a9 | 156 | fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) |
30159e6f | 157 | |
f507c09b | 158 | fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data()); |
30159e6f | 159 | |
160 | fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2); | |
161 | fOutputList->Add(fEvtSel); | |
0b21f686 | 162 | |
163 | fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5); | |
164 | fOutputList->Add(fNClusEt10); | |
30159e6f | 165 | |
cc57d293 | 166 | fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5); |
167 | fOutputList->Add(fRecoPV); | |
168 | ||
30159e6f | 169 | fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100); |
170 | fOutputList->Add(fPVtxZ); | |
c7bb0b43 | 171 | |
f507c09b | 172 | fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5); |
173 | fOutputList->Add(fTrMultDist); | |
174 | ||
caaf99d3 | 175 | fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20); |
176 | fMCDirPhotonPtEtaPhi->Sumw2(); | |
177 | fOutputList->Add(fMCDirPhotonPtEtaPhi); | |
c7bb0b43 | 178 | |
179 | fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100); | |
180 | fDecayPhotonPtMC->Sumw2(); | |
181 | fOutputList->Add(fDecayPhotonPtMC); | |
182 | ||
30159e6f | 183 | 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 | 184 | fOutputList->Add(fCellAbsIdVsAmpl); |
30159e6f | 185 | |
186 | 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); | |
187 | fOutputList->Add(fNClusHighClusE); | |
188 | ||
f9e2362a | 189 | 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); |
190 | fOutputList->Add(fHigherPtConeM02); | |
191 | ||
f507c09b | 192 | fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100); |
193 | fOutputList->Add(fClusEtMcPt); | |
194 | ||
195 | fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7); | |
196 | fOutputList->Add(fClusMcDetaDphi); | |
197 | ||
198 | fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5); | |
199 | fOutputList->Add(fNClusPerPho); | |
200 | ||
d4f449df | 201 | 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 | 202 | fOutputList->Add(fMcPtInConeBG); |
203 | ||
d4f449df | 204 | 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 | 205 | fOutputList->Add(fMcPtInConeSBG); |
206 | ||
d4f449df | 207 | 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 | 208 | fOutputList->Add(fMcPtInConeBGnoUE); |
209 | ||
d4f449df | 210 | 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 | 211 | fOutputList->Add(fMcPtInConeSBGnoUE); |
212 | ||
213 | ||
092ceec8 | 214 | 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); |
215 | fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus); | |
f507c09b | 216 | |
217 | 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=101; | |
218 | Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt}; | |
c1f18270 | 219 | fNDimensions = sizeof(bins)/sizeof(Int_t); |
220 | const Int_t ndims = fNDimensions; | |
f507c09b | 221 | Double_t xmin[] = { 0., 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-1.5}; |
222 | Double_t xmax[] = { 100., 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5}; | |
223 | 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 | 224 | fOutputList->Add(fHnOutput); |
225 | ||
226 | ||
227 | ||
30159e6f | 228 | PostData(1, fOutputList); |
229 | } | |
230 | ||
231 | //________________________________________________________________________ | |
232 | void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) | |
233 | { | |
bd092f0f | 234 | // Main loop, called for each event. |
235 | ||
236 | // event trigger selection | |
30159e6f | 237 | Bool_t isSelected = 0; |
751194e8 | 238 | if(fPeriod.Contains("11a")){ |
239 | if(fTrigBit.Contains("kEMC")) | |
240 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1); | |
241 | else | |
242 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); | |
243 | } | |
244 | else{ | |
245 | if(fTrigBit.Contains("kEMC")) | |
246 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7); | |
247 | else | |
f507c09b | 248 | if(fTrigBit.Contains("kMB")) |
249 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); | |
250 | else | |
251 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); | |
751194e8 | 252 | } |
c7bb0b43 | 253 | if(fIsMc) |
254 | isSelected = kTRUE; | |
ecd47673 | 255 | if(fDebug) |
256 | printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc); | |
30159e6f | 257 | if(!isSelected ) |
258 | return; | |
f3843637 | 259 | if(fIsMc){ |
260 | TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree(); | |
261 | TFile *file = (TFile*)tree->GetCurrentFile(); | |
262 | TString filename = file->GetName(); | |
263 | if(!filename.Contains(fPathStrOpt.Data())) | |
264 | return; | |
265 | } | |
30159e6f | 266 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); |
267 | if (!fESD) { | |
268 | printf("ERROR: fESD not available\n"); | |
269 | return; | |
270 | } | |
271 | ||
272 | fEvtSel->Fill(0); | |
ecd47673 | 273 | if(fDebug) |
274 | printf("fESD is ok\n"); | |
30159e6f | 275 | |
276 | AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex(); | |
d1582697 | 277 | if(!pv) |
30159e6f | 278 | return; |
d1582697 | 279 | if(!pv->GetStatus()) |
280 | fRecoPV->Fill(0); | |
cc57d293 | 281 | else |
282 | fRecoPV->Fill(1); | |
30159e6f | 283 | fPVtxZ->Fill(pv->GetZ()); |
284 | if(TMath::Abs(pv->GetZ())>15) | |
285 | return; | |
ecd47673 | 286 | if(fDebug) |
287 | printf("passed vertex cut\n"); | |
30159e6f | 288 | |
289 | fEvtSel->Fill(1); | |
bd092f0f | 290 | |
06a28959 | 291 | fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks")); |
292 | if(!fTracks){ | |
293 | AliError("Track array in event is NULL!"); | |
294 | return; | |
295 | } | |
30159e6f | 296 | // Track loop to fill a pT spectrum |
3f4073ba | 297 | const Int_t Ntracks = fTracks->GetEntriesFast(); |
298 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { | |
299 | // for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { | |
300 | //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks); | |
301 | AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks)); | |
30159e6f | 302 | if (!track) |
303 | continue; | |
304 | if (fPrTrCuts && fPrTrCuts->IsSelected(track)){ | |
305 | fSelPrimTracks->Add(track); | |
306 | //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi()); | |
307 | } | |
bd092f0f | 308 | } |
309 | ||
30159e6f | 310 | if(!fIsTrain){ |
311 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ | |
312 | if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) | |
313 | break; | |
314 | fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod); | |
315 | } | |
316 | } | |
16a4050e | 317 | AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); |
f507c09b | 318 | fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); |
319 | fTrMultDist->Fill(fTrackMult); | |
16a4050e | 320 | |
30159e6f | 321 | fESD->GetEMCALClusters(fCaloClusters); |
322 | fEMCalCells = fESD->GetEMCALCells(); | |
323 | ||
324 | ||
c7bb0b43 | 325 | |
326 | fMCEvent = MCEvent(); | |
ecd47673 | 327 | if(fMCEvent){ |
328 | if(fDebug) | |
ec9ab86b | 329 | std::cout<<"MCevent exists\n"; |
c7bb0b43 | 330 | fStack = (AliStack*)fMCEvent->Stack(); |
ecd47673 | 331 | } |
332 | else{ | |
333 | if(fDebug) | |
ec9ab86b | 334 | std::cout<<"ERROR: NO MC EVENT!!!!!!\n"; |
ecd47673 | 335 | } |
f507c09b | 336 | FollowGamma(); |
337 | if(fDebug) | |
338 | printf("passed calling of FollowGamma\n"); | |
339 | FillClusHists(); | |
340 | if(fDebug) | |
341 | printf("passed calling of FillClusHists\n"); | |
63e40cb8 | 342 | FillMcHists(); |
ecd47673 | 343 | if(fDebug) |
344 | printf("passed calling of FillMcHists\n"); | |
345 | ||
f507c09b | 346 | fCaloClusters->Clear(); |
347 | fSelPrimTracks->Clear(); | |
348 | fNClusForDirPho = 0; | |
c7bb0b43 | 349 | |
30159e6f | 350 | PostData(1, fOutputList); |
351 | } | |
bd092f0f | 352 | |
30159e6f | 353 | //________________________________________________________________________ |
354 | void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() | |
355 | { | |
bd092f0f | 356 | // Fill cluster histograms. |
357 | ||
30159e6f | 358 | if(!fCaloClusters) |
359 | return; | |
360 | const Int_t nclus = fCaloClusters->GetEntries(); | |
361 | if(nclus==0) | |
362 | return; | |
ecd47673 | 363 | if(fDebug) |
364 | printf("Inside FillClusHists and there are %d clusters\n",nclus); | |
e681d9ce | 365 | Double_t maxE = 0; |
0b21f686 | 366 | Int_t nclus10 = 0; |
f507c09b | 367 | Double_t ptmc=-1; |
30159e6f | 368 | for(Int_t ic=0;ic<nclus;ic++){ |
369 | maxE=0; | |
370 | AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic)); | |
371 | if(!c) | |
372 | continue; | |
373 | if(!c->IsEMCAL()) | |
374 | continue; | |
965c985f | 375 | if(c->E()<fECut) |
376 | continue; | |
30159e6f | 377 | Short_t id; |
378 | Double_t Emax = GetMaxCellEnergy( c, id); | |
379 | Double_t Ecross = GetCrossEnergy( c, id); | |
380 | if((1-Ecross/Emax)>fExoticCut) | |
381 | continue; | |
382 | Float_t clsPos[3] = {0,0,0}; | |
383 | c->GetPosition(clsPos); | |
384 | TVector3 clsVec(clsPos); | |
3dc2825e | 385 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); |
f507c09b | 386 | if(fDebug) |
387 | printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E()); | |
0b21f686 | 388 | if(Et>10) |
389 | nclus10++; | |
30159e6f | 390 | Float_t ceiso, cephiband, cecore; |
391 | Float_t triso, trphiband, trcore; | |
392 | Float_t alliso, allphiband, allcore; | |
393 | Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR; | |
394 | Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04); | |
f224025f | 395 | GetCeIso(clsVec, id, ceiso, cephiband, cecore); |
30159e6f | 396 | GetTrIso(clsVec, triso, trphiband, trcore); |
f9e2362a | 397 | Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz()); |
398 | if(Et>10 && Et<15 && dr>0.025){ | |
399 | fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02()); | |
400 | if(fDebug) | |
401 | printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02()); | |
402 | } | |
30159e6f | 403 | alliso = ceiso + triso; |
404 | allphiband = cephiband + trphiband; | |
405 | allcore = cecore + trcore; | |
406 | Float_t ceisoue = cephiband/phibandArea*netConeArea; | |
407 | Float_t trisoue = trphiband/phibandArea*netConeArea; | |
408 | Float_t allisoue = allphiband/phibandArea*netConeArea; | |
f912f9a9 | 409 | Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR); |
410 | if(fDebug && Et>10) | |
411 | printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et); | |
412 | if(c->GetM02()>0.5 && c->GetM02()<2.0){ | |
413 | fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum); | |
414 | fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum); | |
415 | } | |
416 | if(c->GetM02()>0.1 && c->GetM02()<0.3){ | |
417 | fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum); | |
418 | fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum); | |
419 | } | |
965c985f | 420 | const Int_t ndims = fNDimensions; |
421 | Double_t outputValues[ndims]; | |
f507c09b | 422 | ptmc = GetClusSource(c); |
965c985f | 423 | outputValues[0] = Et; |
424 | outputValues[1] = c->GetM02(); | |
f912f9a9 | 425 | outputValues[2] = ceiso-Et/*cecore*/-ceisoue; |
965c985f | 426 | outputValues[3] = triso-trisoue; |
f912f9a9 | 427 | outputValues[4] = alliso-Et/*cecore*/-allisoue; |
302b398c | 428 | outputValues[5] = ceiso-Et; |
429 | outputValues[6] = alliso-Et; | |
717bb45b | 430 | outputValues[7] = c->GetTrackDx(); |
431 | outputValues[8] = c->GetTrackDz(); | |
432 | outputValues[9] = clsVec.Eta(); | |
433 | outputValues[10] = clsVec.Phi(); | |
16a4050e | 434 | outputValues[11] = fEMCalCells->GetCellTime(id); |
435 | outputValues[12] = fTrackMult; | |
f507c09b | 436 | outputValues[13] = ptmc; |
965c985f | 437 | fHnOutput->Fill(outputValues); |
30159e6f | 438 | if(c->E()>maxE) |
439 | maxE = c->E(); | |
440 | } | |
441 | fNClusHighClusE->Fill(maxE,nclus); | |
0b21f686 | 442 | fNClusEt10->Fill(nclus10); |
f507c09b | 443 | fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho); |
30159e6f | 444 | } |
bd092f0f | 445 | |
30159e6f | 446 | //________________________________________________________________________ |
f224025f | 447 | void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core) |
30159e6f | 448 | { |
bd092f0f | 449 | // Get cell isolation. |
450 | ||
30159e6f | 451 | if(!fEMCalCells) |
452 | return; | |
453 | const Int_t ncells = fEMCalCells->GetNumberOfCells(); | |
454 | Float_t totiso=0; | |
455 | Float_t totband=0; | |
456 | Float_t totcore=0; | |
457 | Float_t etacl = vec.Eta(); | |
458 | Float_t phicl = vec.Phi(); | |
459 | Float_t thetacl = vec.Theta(); | |
f224025f | 460 | Float_t maxtcl = fEMCalCells->GetCellTime(maxid); |
30159e6f | 461 | if(phicl<0) |
462 | phicl+=TMath::TwoPi(); | |
463 | Int_t absid = -1; | |
464 | Float_t eta=-1, phi=-1; | |
465 | for(int icell=0;icell<ncells;icell++){ | |
466 | absid = TMath::Abs(fEMCalCells->GetCellNumber(icell)); | |
c3f33463 | 467 | Float_t celltime = fEMCalCells->GetCellTime(absid); |
f224025f | 468 | //if(TMath::Abs(celltime)>2e-8 && (!fIsMc)) |
469 | if(TMath::Abs(celltime-maxtcl)>2e-8 ) | |
c3f33463 | 470 | continue; |
30159e6f | 471 | if(!fGeom) |
472 | return; | |
473 | fGeom->EtaPhiFromIndex(absid,eta,phi); | |
474 | Float_t dphi = TMath::Abs(phi-phicl); | |
475 | Float_t deta = TMath::Abs(eta-etacl); | |
476 | Float_t R = TMath::Sqrt(deta*deta + dphi*dphi); | |
477 | Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl); | |
478 | if(R<fIsoConeR){ | |
479 | totiso += etcell; | |
480 | if(R<0.04) | |
481 | totcore += etcell; | |
482 | } | |
483 | else{ | |
484 | if(dphi>fIsoConeR) | |
485 | continue; | |
486 | if(deta<fIsoConeR) | |
487 | continue; | |
488 | totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl); | |
489 | } | |
490 | } | |
491 | iso = totiso; | |
492 | phiband = totband; | |
493 | core = totcore; | |
494 | } | |
bd092f0f | 495 | |
30159e6f | 496 | //________________________________________________________________________ |
497 | void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core) | |
498 | { | |
bd092f0f | 499 | // Get track isolation. |
500 | ||
30159e6f | 501 | if(!fSelPrimTracks) |
502 | return; | |
f9e2362a | 503 | fHigherPtCone = 0; |
30159e6f | 504 | const Int_t ntracks = fSelPrimTracks->GetEntries(); |
505 | Double_t totiso=0; | |
506 | Double_t totband=0; | |
507 | Double_t totcore=0; | |
508 | Double_t etacl = vec.Eta(); | |
509 | Double_t phicl = vec.Phi(); | |
510 | if(phicl<0) | |
511 | phicl+=TMath::TwoPi(); | |
512 | for(int itrack=0;itrack<ntracks;itrack++){ | |
3f4073ba | 513 | AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack)); |
30159e6f | 514 | if(!track) |
515 | continue; | |
516 | Double_t dphi = TMath::Abs(track->Phi()-phicl); | |
517 | Double_t deta = TMath::Abs(track->Eta()-etacl); | |
518 | Double_t R = TMath::Sqrt(deta*deta + dphi*dphi); | |
519 | Double_t pt = track->Pt(); | |
f9e2362a | 520 | if(pt>fHigherPtCone) |
521 | fHigherPtCone = pt; | |
30159e6f | 522 | if(R<fIsoConeR){ |
523 | totiso += track->Pt(); | |
524 | if(R<0.04) | |
525 | totcore += pt; | |
526 | } | |
527 | else{ | |
528 | if(dphi>fIsoConeR) | |
529 | continue; | |
530 | if(deta<fIsoConeR) | |
531 | continue; | |
532 | totband += track->Pt(); | |
533 | } | |
534 | } | |
535 | iso = totiso; | |
536 | phiband = totband; | |
537 | core = totcore; | |
538 | } | |
bd092f0f | 539 | |
30159e6f | 540 | //________________________________________________________________________ |
541 | Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) | |
542 | { | |
543 | // Calculate the energy of cross cells around the leading cell. | |
544 | ||
545 | AliVCaloCells *cells = 0; | |
546 | cells = fESD->GetEMCALCells(); | |
547 | if (!cells) | |
548 | return 0; | |
549 | ||
30159e6f | 550 | if (!fGeom) |
551 | return 0; | |
552 | ||
553 | Int_t iSupMod = -1; | |
554 | Int_t iTower = -1; | |
555 | Int_t iIphi = -1; | |
556 | Int_t iIeta = -1; | |
557 | Int_t iphi = -1; | |
558 | Int_t ieta = -1; | |
559 | Int_t iphis = -1; | |
560 | Int_t ietas = -1; | |
561 | ||
562 | Double_t crossEnergy = 0; | |
563 | ||
564 | fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta); | |
565 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas); | |
566 | ||
567 | Int_t ncells = cluster->GetNCells(); | |
568 | for (Int_t i=0; i<ncells; i++) { | |
569 | Int_t cellAbsId = cluster->GetCellAbsId(i); | |
570 | fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); | |
571 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); | |
572 | Int_t aphidiff = TMath::Abs(iphi-iphis); | |
573 | if (aphidiff>1) | |
574 | continue; | |
575 | Int_t aetadiff = TMath::Abs(ieta-ietas); | |
576 | if (aetadiff>1) | |
577 | continue; | |
578 | if ( (aphidiff==1 && aetadiff==0) || | |
579 | (aphidiff==0 && aetadiff==1) ) { | |
580 | crossEnergy += cells->GetCellAmplitude(cellAbsId); | |
581 | } | |
582 | } | |
583 | ||
584 | return crossEnergy; | |
585 | } | |
586 | ||
30159e6f | 587 | //________________________________________________________________________ |
588 | Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const | |
589 | { | |
590 | // Get maximum energy of attached cell. | |
591 | ||
592 | id = -1; | |
593 | ||
594 | AliVCaloCells *cells = 0; | |
595 | cells = fESD->GetEMCALCells(); | |
596 | if (!cells) | |
597 | return 0; | |
598 | ||
599 | Double_t maxe = 0; | |
600 | Int_t ncells = cluster->GetNCells(); | |
601 | for (Int_t i=0; i<ncells; i++) { | |
602 | Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); | |
603 | if (e>maxe) { | |
604 | maxe = e; | |
605 | id = cluster->GetCellAbsId(i); | |
606 | } | |
607 | } | |
608 | return maxe; | |
609 | } | |
bd092f0f | 610 | |
c7bb0b43 | 611 | //________________________________________________________________________ |
612 | void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists() | |
613 | { | |
614 | if(!fStack) | |
615 | return; | |
616 | Int_t nTracks = fStack->GetNtrack(); | |
ecd47673 | 617 | if(fDebug) |
618 | printf("Inside FillMcHists and there are %d mcparts\n",nTracks); | |
c7bb0b43 | 619 | for(Int_t iTrack=0;iTrack<nTracks;iTrack++){ |
620 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); | |
621 | if(!mcp) | |
622 | continue; | |
623 | Int_t pdg = mcp->GetPdgCode(); | |
624 | if(pdg!=22) | |
625 | continue; | |
f2acdbe9 | 626 | if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2) |
627 | continue; | |
c7bb0b43 | 628 | Int_t imom = mcp->GetMother(0); |
629 | if(imom<0 || imom>nTracks) | |
630 | continue; | |
631 | TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom)); | |
632 | if(!mcmom) | |
633 | continue; | |
634 | Int_t pdgMom = mcmom->GetPdgCode(); | |
ec9ab86b | 635 | if((imom==6 || imom==7) && pdgMom==22) { |
caaf99d3 | 636 | fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
f507c09b | 637 | if(fNClusForDirPho==0) |
092ceec8 | 638 | fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
40c4f1bf | 639 | if(fDebug){ |
f2acdbe9 | 640 | 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 | 641 | 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()); |
642 | } | |
f2acdbe9 | 643 | } |
c7bb0b43 | 644 | else{ |
645 | if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000) | |
646 | fDecayPhotonPtMC->Fill(mcp->Pt()); | |
647 | } | |
648 | } | |
649 | } | |
30159e6f | 650 | //________________________________________________________________________ |
f507c09b | 651 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) |
652 | { | |
653 | if(!c) | |
654 | return -1; | |
655 | if(!fStack) | |
656 | return -1; | |
657 | Int_t nmcp = fStack->GetNtrack(); | |
658 | Int_t clabel = c->GetLabel(); | |
659 | if(fDebug && fMcIdFamily.Contains(Form("%d",clabel))) | |
660 | printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel); | |
661 | if(!fMcIdFamily.Contains(Form("%d",clabel))) | |
662 | return -1; | |
663 | fNClusForDirPho++; | |
664 | TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1); | |
665 | Int_t partonpos = partonposstr.Atoi(); | |
666 | if(fDebug) | |
667 | printf("\t^^^^ parton position = %d ^^^^\n",partonpos); | |
668 | if(clabel<0 || clabel>nmcp) | |
669 | return 0; | |
670 | Float_t clsPos[3] = {0,0,0}; | |
671 | c->GetPosition(clsPos); | |
672 | TVector3 clsVec(clsPos); | |
673 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); | |
674 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos)); | |
675 | if(!mcp) | |
676 | return -1; | |
677 | if(fDebug){ | |
678 | 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); | |
679 | } | |
680 | Int_t lab1 = mcp->GetFirstDaughter(); | |
681 | if(lab1<0 || lab1>nmcp) | |
682 | return -1; | |
683 | TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1)); | |
684 | if(!mcd) | |
685 | return -1; | |
686 | if(fDebug) | |
687 | 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); | |
688 | if(mcd->GetPdgCode()==22){ | |
689 | fClusEtMcPt->Fill(mcd->Pt(), Et); | |
690 | fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi()); | |
691 | } | |
692 | else{ | |
693 | printf("Warning: daughter of photon parton is not a photon\n"); | |
694 | return -1; | |
695 | } | |
696 | return fDirPhoPt; | |
697 | } | |
698 | //________________________________________________________________________ | |
699 | void AliAnalysisTaskEMCALIsoPhoton::FollowGamma() | |
700 | { | |
701 | if(!fStack) | |
702 | return; | |
703 | Int_t selfid = 6; | |
704 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid)); | |
705 | if(!mcp) | |
706 | return; | |
707 | if(mcp->GetPdgCode()!=22){ | |
708 | mcp = static_cast<TParticle*>(fStack->Particle(++selfid)); | |
709 | if(!mcp) | |
710 | return; | |
711 | } | |
712 | Int_t daug0f = mcp->GetFirstDaughter(); | |
713 | Int_t daug0l = mcp->GetLastDaughter(); | |
714 | Int_t nd0 = daug0l - daug0f; | |
715 | if(fDebug) | |
716 | 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); | |
717 | fMcIdFamily = Form("%d,",selfid); | |
718 | GetDaughtersInfo(daug0f,daug0l, selfid,""); | |
719 | if(fDebug){ | |
720 | printf("\t---- end of the prompt gamma's family tree ----\n\n"); | |
721 | printf("Family id string = %s,\n\n",fMcIdFamily.Data()); | |
722 | } | |
723 | TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f)); | |
724 | if(!md) | |
725 | return; | |
726 | fDirPhoPt = md->Pt(); | |
727 | } | |
728 | //________________________________________________________________________ | |
22ad7981 | 729 | void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind) |
f507c09b | 730 | { |
731 | int nmcp = fStack->GetNtrack(); | |
732 | if(firstd<0 || firstd>nmcp) | |
733 | return; | |
734 | if(lastd<0 || lastd>nmcp) | |
735 | return; | |
736 | if(firstd>lastd){ | |
737 | printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd); | |
738 | return; | |
739 | } | |
740 | TString indenter = Form("\t%s",inputind); | |
741 | TParticle *dp = 0x0; | |
742 | if(fDebug) | |
743 | printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid); | |
744 | for(int id=firstd; id<lastd+1; id++){ | |
745 | dp = static_cast<TParticle*>(fStack->Particle(id)); | |
746 | if(!dp) | |
747 | continue; | |
748 | Int_t dfd = dp->GetFirstDaughter(); | |
749 | Int_t dld = dp->GetLastDaughter(); | |
750 | Int_t dnd = dld - dfd + 1; | |
751 | Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy()); | |
752 | if(fDebug) | |
753 | 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); | |
754 | fMcIdFamily += Form("%d,",id); | |
755 | GetDaughtersInfo(dfd,dld,id,indenter.Data()); | |
756 | } | |
757 | } | |
f912f9a9 | 758 | |
759 | //________________________________________________________________________ | |
760 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R) | |
761 | { | |
762 | if(!fStack) | |
763 | return 0; | |
764 | if(fDebug) | |
765 | printf("Inside GetMcPtSumInCone!!\n"); | |
766 | Int_t nTracks = fStack->GetNtrack(); | |
767 | Float_t ptsum = 0; | |
768 | for(Int_t iTrack=0;iTrack<nTracks;iTrack++){ | |
769 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); | |
770 | if(!mcp) | |
771 | continue; | |
772 | if(mcp->Rho()>2.5) | |
773 | continue; | |
774 | else { | |
775 | if(fDebug) | |
776 | printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz()); | |
777 | } | |
778 | Float_t dphi = mcp->Phi() - phiclus; | |
779 | Float_t deta = mcp->Eta() - etaclus; | |
780 | if(fDebug) | |
781 | printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta); | |
782 | if(deta>10) | |
783 | continue; | |
784 | Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta); | |
785 | if(dR>R) | |
786 | continue; | |
787 | ptsum += mcp->Pt(); | |
788 | } | |
789 | return ptsum; | |
790 | } | |
f507c09b | 791 | //________________________________________________________________________ |
30159e6f | 792 | void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) |
793 | { | |
bd092f0f | 794 | // Called once at the end of the query. |
30159e6f | 795 | } |