48d874ff |
1 | // $Id$ |
542ac102 |
2 | // |
3 | // Emcal jet finder task. |
4 | // |
b65fac7a |
5 | // Authors: C.Loizides, S.Aiola |
48d874ff |
6 | |
c64cb1f6 |
7 | #include <vector> |
04905231 |
8 | #include "AliEmcalJetTask.h" |
9 | |
6ea168d0 |
10 | #include <TChain.h> |
48d874ff |
11 | #include <TClonesArray.h> |
6ea168d0 |
12 | #include <TList.h> |
13 | #include <TLorentzVector.h> |
b65fac7a |
14 | |
48d874ff |
15 | #include "AliAnalysisManager.h" |
6ea168d0 |
16 | #include "AliCentrality.h" |
f5f3c8e9 |
17 | #include "AliEMCALGeometry.h" |
04905231 |
18 | #include "AliESDEvent.h" |
f5f3c8e9 |
19 | #include "AliEmcalJet.h" |
04905231 |
20 | #include "AliEmcalParticle.h" |
63206072 |
21 | #include "AliFJWrapper.h" |
04905231 |
22 | #include "AliMCEvent.h" |
35789a2d |
23 | #include "AliVCluster.h" |
b65fac7a |
24 | #include "AliVEvent.h" |
04905231 |
25 | #include "AliVParticle.h" |
48d874ff |
26 | |
914d486c |
27 | ClassImp(AliEmcalJetTask) |
48d874ff |
28 | |
48d874ff |
29 | //________________________________________________________________________ |
914d486c |
30 | AliEmcalJetTask::AliEmcalJetTask() : |
31 | AliAnalysisTaskSE("AliEmcalJetTask"), |
6ea168d0 |
32 | fTracksName("Tracks"), |
33 | fCaloName("CaloClusters"), |
34 | fJetsName("Jets"), |
35 | fAlgo(1), |
36 | fRadius(0.4), |
37 | fType(0), |
914d486c |
38 | fMinJetTrackPt(0.15), |
39 | fMinJetClusPt(0.15), |
04905231 |
40 | fPhiMin(0), |
41 | fPhiMax(TMath::TwoPi()), |
42 | fEtaMin(-1), |
43 | fEtaMax(1), |
f5f3c8e9 |
44 | fMinJetArea(0.01), |
45 | fMinJetPt(1.0), |
04905231 |
46 | fIsInit(0), |
47 | fIsMcPart(0), |
48 | fIsEmcPart(0), |
b65fac7a |
49 | fJets(0), |
04905231 |
50 | fEvent(0), |
51 | fTracks(0), |
52 | fClus(0) |
6ea168d0 |
53 | { |
54 | // Default constructor. |
6ea168d0 |
55 | } |
56 | |
57 | //________________________________________________________________________ |
914d486c |
58 | AliEmcalJetTask::AliEmcalJetTask(const char *name) : |
63206072 |
59 | AliAnalysisTaskSE(name), |
7efbea04 |
60 | fTracksName("Tracks"), |
48d874ff |
61 | fCaloName("CaloClusters"), |
7efbea04 |
62 | fJetsName("Jets"), |
48d874ff |
63 | fAlgo(1), |
64 | fRadius(0.4), |
65 | fType(0), |
6ea168d0 |
66 | fMinJetTrackPt(0.15), |
914d486c |
67 | fMinJetClusPt(0.15), |
04905231 |
68 | fPhiMin(0), |
69 | fPhiMax(TMath::TwoPi()), |
70 | fEtaMin(-1), |
71 | fEtaMax(1), |
f5f3c8e9 |
72 | fMinJetArea(0.01), |
73 | fMinJetPt(1.0), |
04905231 |
74 | fIsInit(0), |
75 | fIsMcPart(0), |
76 | fIsEmcPart(0), |
b65fac7a |
77 | fJets(0), |
04905231 |
78 | fEvent(0), |
79 | fTracks(0), |
80 | fClus(0) |
48d874ff |
81 | { |
82 | // Standard constructor. |
7efbea04 |
83 | |
7efbea04 |
84 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex."; |
48d874ff |
85 | } |
86 | |
87 | //________________________________________________________________________ |
914d486c |
88 | AliEmcalJetTask::~AliEmcalJetTask() |
48d874ff |
89 | { |
90 | // Destructor |
91 | } |
92 | |
93 | //________________________________________________________________________ |
914d486c |
94 | void AliEmcalJetTask::UserCreateOutputObjects() |
48d874ff |
95 | { |
96 | // Create user objects. |
97 | |
914d486c |
98 | fJets = new TClonesArray("AliEmcalJet"); |
48d874ff |
99 | fJets->SetName(fJetsName); |
48d874ff |
100 | } |
101 | |
102 | //________________________________________________________________________ |
914d486c |
103 | void AliEmcalJetTask::UserExec(Option_t *) |
48d874ff |
104 | { |
105 | // Main loop, called for each event. |
106 | |
04905231 |
107 | if (!fIsInit) { |
108 | if (!DoInit()) |
109 | return; |
110 | fIsInit = kTRUE; |
1772fe65 |
111 | } |
b65fac7a |
112 | |
63206072 |
113 | // delete jet output |
114 | fJets->Delete(); |
115 | |
04905231 |
116 | FindJets(); |
48d874ff |
117 | } |
118 | |
119 | //________________________________________________________________________ |
914d486c |
120 | void AliEmcalJetTask::Terminate(Option_t *) |
48d874ff |
121 | { |
122 | // Called once at the end of the analysis. |
48d874ff |
123 | } |
124 | |
125 | //________________________________________________________________________ |
04905231 |
126 | void AliEmcalJetTask::FindJets() |
48d874ff |
127 | { |
128 | // Find jets. |
129 | |
04905231 |
130 | if (!fTracks && !fClus) |
d03084cd |
131 | return; |
132 | |
48d874ff |
133 | TString name("kt"); |
134 | fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm); |
04905231 |
135 | if (fAlgo>=1) { |
48d874ff |
136 | name = "antikt"; |
137 | jalgo = fastjet::antikt_algorithm; |
138 | } |
139 | |
04905231 |
140 | // setup fj wrapper |
48d874ff |
141 | AliFJWrapper fjw(name, name); |
f5f3c8e9 |
142 | fjw.SetAreaType(fastjet::active_area_explicit_ghosts); |
04905231 |
143 | fjw.SetR(fRadius); |
6ea168d0 |
144 | fjw.SetAlgorithm(jalgo); |
04905231 |
145 | fjw.SetMaxRap(1); |
48d874ff |
146 | fjw.Clear(); |
147 | |
04905231 |
148 | // get primary vertex |
149 | Double_t vertex[3] = {0, 0, 0}; |
150 | fEvent->GetPrimaryVertex()->GetXYZ(vertex); |
151 | |
152 | if (fTracks) { |
153 | const Int_t Ntracks = fTracks->GetEntries(); |
48d874ff |
154 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { |
04905231 |
155 | AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks)); |
48d874ff |
156 | if (!t) |
157 | continue; |
1772fe65 |
158 | if (t->Pt() < fMinJetTrackPt) |
914d486c |
159 | continue; |
04905231 |
160 | Double_t eta = t->Eta(); |
161 | Double_t phi = t->Phi(); |
162 | if ((eta<fEtaMin) && (eta>fEtaMax) && |
163 | (phi<fPhiMin) && (phi<fPhiMax)) |
164 | continue; |
a3c8c8c8 |
165 | // offset of 100 for consistency with cluster ids |
166 | fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100); |
48d874ff |
167 | } |
168 | } |
6ea168d0 |
169 | |
04905231 |
170 | if (fClus) { |
171 | if (fIsEmcPart) { |
172 | const Int_t Nclus = fClus->GetEntries(); |
173 | for (Int_t iClus = 0; iClus < Nclus; ++iClus) { |
174 | AliVCluster *c = 0; |
175 | Double_t cEta=0,cPhi=0,cPt=0; |
176 | Double_t cPx=0,cPy=0,cPz=0; |
177 | if (fIsEmcPart) { |
178 | AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus)); |
179 | if (!ep) |
180 | continue; |
181 | c = ep->GetCluster(); |
182 | if (!c) |
183 | continue; |
184 | cEta = ep->Eta(); |
185 | cPhi = ep->Phi(); |
186 | cPt = ep->Pt(); |
187 | cPx = ep->Px(); |
188 | cPy = ep->Py(); |
189 | cPz = ep->Pz(); |
190 | } else { |
191 | c = static_cast<AliVCluster*>(fClus->At(iClus)); |
192 | if (!c) |
193 | continue; |
194 | TLorentzVector nP; |
195 | c->GetMomentum(nP, vertex); |
196 | cEta = nP.Eta(); |
197 | cPhi = nP.Phi(); |
198 | cPt = nP.Pt(); |
199 | cPx = nP.Px(); |
200 | cPy = nP.Py(); |
201 | cPz = nP.Pz(); |
202 | } |
203 | if (!c->IsEMCAL()) |
204 | continue; |
205 | if (cPt < fMinJetClusPt) |
206 | continue; |
207 | if ((cEta<fEtaMin) && (cEta>fEtaMax) && |
208 | (cPhi<fPhiMin) && (cPhi<fPhiMax)) |
209 | continue; |
210 | // offset of 100 to skip ghost particles uid = -1 |
211 | fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100); |
212 | } |
213 | } else { /*CaloClusters given as input*/ |
214 | const Int_t Nclus = fClus->GetEntries(); |
215 | for (Int_t iClus = 0; iClus < Nclus; ++iClus) { |
216 | AliVCluster *c = static_cast<AliVCluster*>(fClus->At(iClus)); |
217 | if (!c) |
218 | continue; |
219 | if (!c->IsEMCAL()) |
220 | continue; |
221 | TLorentzVector nPart; |
222 | c->GetMomentum(nPart, vertex); |
223 | if (nPart.Pt() < fMinJetClusPt) |
224 | continue; |
225 | // offset of 100 to skip ghost particles uid = -1 |
226 | fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100); |
227 | } |
48d874ff |
228 | } |
229 | } |
04905231 |
230 | |
7efbea04 |
231 | // run jet finder |
48d874ff |
232 | fjw.Run(); |
f5f3c8e9 |
233 | |
234 | // get geometry |
235 | AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); |
236 | if (!geom) { |
02c7e943 |
237 | AliFatal(Form("%s: Can not create geometry", GetName())); |
f5f3c8e9 |
238 | return; |
239 | } |
04905231 |
240 | |
241 | // loop over fastjet jets |
48d874ff |
242 | std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets(); |
b65fac7a |
243 | for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) { |
f5f3c8e9 |
244 | if (jets_incl[ij].perp()<fMinJetPt) |
245 | continue; |
246 | if (fjw.GetJetArea(ij)<fMinJetArea) |
48d874ff |
247 | continue; |
04905231 |
248 | |
914d486c |
249 | AliEmcalJet *jet = new ((*fJets)[jetCount]) |
250 | AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m()); |
04905231 |
251 | |
252 | // loop over constituents |
c64cb1f6 |
253 | std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij)); |
04905231 |
254 | jet->SetNumberOfTracks(constituents.size()); |
255 | jet->SetNumberOfClusters(constituents.size()); |
256 | |
b65fac7a |
257 | Int_t nt = 0; |
258 | Int_t nc = 0; |
f5f3c8e9 |
259 | Double_t neutralE = 0; |
111318e8 |
260 | Double_t maxCh = 0; |
261 | Double_t maxNe = 0; |
f5f3c8e9 |
262 | Int_t gall = 0; |
263 | Int_t gemc = 0; |
04905231 |
264 | Int_t cemc = 0; |
111318e8 |
265 | Int_t ncharged = 0; |
266 | Int_t nneutral = 0; |
02c7e943 |
267 | Double_t mcpt = 0; |
04905231 |
268 | Double_t emcpt = 0; |
b65fac7a |
269 | |
b65fac7a |
270 | for(UInt_t ic = 0; ic < constituents.size(); ++ic) { |
48d874ff |
271 | Int_t uid = constituents[ic].user_index(); |
b65fac7a |
272 | |
04905231 |
273 | if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle |
f5f3c8e9 |
274 | ++gall; |
04905231 |
275 | Double_t gphi = constituents[ic].phi(); |
276 | if (gphi<0) |
277 | gphi += TMath::TwoPi(); |
278 | gphi *= TMath::RadToDeg(); |
f5f3c8e9 |
279 | Double_t geta = constituents[ic].eta(); |
b65fac7a |
280 | if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) && |
281 | (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax())) |
04905231 |
282 | ++gemc; |
283 | } else if ((uid > 0) && fTracks) { // track constituent |
111318e8 |
284 | Int_t tid = uid - 100; |
04905231 |
285 | AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid)); |
286 | if (!t) |
287 | continue; |
288 | Double_t cEta = t->Eta(); |
289 | Double_t cPhi = t->Phi(); |
290 | Double_t cPt = t->Pt(); |
291 | Double_t cP = t->P(); |
292 | if (t->Charge() == 0) { |
293 | neutralE += cP; |
294 | ++nneutral; |
295 | if (cPt > maxNe) |
296 | maxNe = cPt; |
297 | } else { |
298 | ++ncharged; |
299 | if (cPt > maxCh) |
300 | maxCh = cPt; |
301 | } |
302 | |
303 | if (fIsMcPart || t->GetLabel() == 100) // check if MC particle |
304 | mcpt += cPt; |
305 | |
306 | if (cPhi<0) |
307 | cPhi += TMath::TwoPi(); |
308 | cPhi *= TMath::RadToDeg(); |
309 | if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) && |
310 | (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) { |
311 | emcpt += cPt; |
312 | ++cemc; |
d03084cd |
313 | } |
04905231 |
314 | |
315 | jet->AddTrackAt(tid, nt); |
316 | ++nt; |
317 | } else if (fClus) { // cluster constituent |
111318e8 |
318 | Int_t cid = -uid - 100; |
04905231 |
319 | AliVCluster *c = 0; |
320 | Double_t cEta=0,cPhi=0,cPt=0,cP=0; |
321 | if (fIsEmcPart) { |
322 | AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid)); |
323 | if (!ep) |
324 | continue; |
325 | c = ep->GetCluster(); |
326 | if (!c) |
327 | continue; |
328 | cEta = ep->Eta(); |
329 | cPhi = ep->Phi(); |
330 | cPt = ep->Pt(); |
331 | cP = ep->P(); |
332 | } else { |
333 | c = static_cast<AliVCluster*>(fClus->At(cid)); |
334 | if (!c) |
335 | continue; |
336 | TLorentzVector nP; |
337 | c->GetMomentum(nP, vertex); |
338 | cEta = nP.Eta(); |
339 | cPhi = nP.Phi(); |
340 | cPt = nP.Pt(); |
341 | cP = nP.P(); |
342 | } |
343 | |
344 | neutralE += cP; |
345 | if (cPt > maxNe) |
346 | maxNe = cPt; |
347 | |
348 | if (c->Chi2() == 100) // MC particle |
349 | mcpt += cPt; |
350 | |
351 | if (cPhi<0) |
352 | cPhi += TMath::TwoPi(); |
353 | cPhi *= TMath::RadToDeg(); |
354 | if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) && |
355 | (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) { |
356 | emcpt += cPt; |
357 | ++cemc; |
358 | } |
359 | |
360 | jet->AddClusterAt(cid, nc); |
361 | ++nc; |
362 | ++nneutral; |
02c7e943 |
363 | } else { |
364 | AliError(Form("%s: No logical way to end up here.", GetName())); |
365 | continue; |
48d874ff |
366 | } |
04905231 |
367 | } /* end of constituent loop */ |
368 | |
35789a2d |
369 | jet->SetNumberOfTracks(nt); |
370 | jet->SetNumberOfClusters(nc); |
f5f3c8e9 |
371 | jet->SortConstituents(); |
111318e8 |
372 | jet->SetMaxNeutralPt(maxNe); |
373 | jet->SetMaxChargedPt(maxCh); |
b65fac7a |
374 | jet->SetNEF(neutralE / jet->E()); |
f5f3c8e9 |
375 | jet->SetArea(fjw.GetJetArea(ij)); |
111318e8 |
376 | jet->SetNumberOfCharged(ncharged); |
377 | jet->SetNumberOfNeutrals(nneutral); |
02c7e943 |
378 | jet->SetMCPt(mcpt); |
04905231 |
379 | jet->SetNEmc(cemc); |
380 | jet->SetPtEmc(emcpt); |
381 | |
b65fac7a |
382 | if (gall > 0) |
383 | jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall); |
f5f3c8e9 |
384 | else |
385 | jet->SetAreaEmc(-1); |
b65fac7a |
386 | if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && |
1772fe65 |
387 | (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) && |
388 | (jet->Eta() > geom->GetArm1EtaMin()) && |
389 | (jet->Eta() < geom->GetArm1EtaMax())) |
f5f3c8e9 |
390 | jet->SetAxisInEmcal(kTRUE); |
48d874ff |
391 | jetCount++; |
392 | } |
393 | } |
04905231 |
394 | |
395 | //________________________________________________________________________ |
396 | Bool_t AliEmcalJetTask::DoInit() |
397 | { |
398 | // Init. Return true if successful. |
399 | |
400 | // get input collections |
401 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); |
402 | |
403 | // get the event |
404 | fEvent = InputEvent(); |
405 | if (!fEvent) { |
406 | AliError(Form("%s: Could not retrieve event! Returning", GetName())); |
407 | return 0; |
408 | } |
409 | |
410 | // add jets to event if not yet there |
411 | if (!(fEvent->FindListObject(fJetsName))) |
412 | fEvent->AddObject(fJets); |
413 | else { |
414 | AliError(Form("%s: Object with name %s already in event! Returning", fJetsName.Data(), GetName())); |
415 | return 0; |
416 | } |
417 | |
418 | if ((fType == 0) || (fType == 1)) { |
419 | if (fTracksName == "Tracks") |
420 | am->LoadBranch("Tracks"); |
421 | if (!fTracks && !fTracksName.IsNull()) { |
422 | fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName)); |
423 | if (!fTracks) { |
424 | AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data())); |
425 | return 0; |
426 | } |
427 | } |
428 | TString objname(fTracks->GetClass()->GetName()); |
429 | TClass cls(objname); |
430 | if (cls.InheritsFrom("AliMCParticle")) |
431 | fIsMcPart = 1; |
432 | } |
433 | |
434 | if ((fType == 0) || (fType == 2)) { |
435 | if (fCaloName == "CaloClusters") |
436 | am->LoadBranch("CaloClusters"); |
437 | if (!fClus && !fCaloName.IsNull()) { |
438 | fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName)); |
439 | if (!fClus) { |
440 | AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data())); |
441 | return 0; |
442 | } |
443 | } |
444 | TString objname(fClus->GetClass()->GetName()); |
445 | TClass cls(objname); |
446 | if (cls.InheritsFrom("AliEmcalParticle")) |
447 | fIsEmcPart = 1; |
448 | } |
449 | |
450 | return 1; |
451 | } |