]>
Commit | Line | Data |
---|---|---|
1 | // $Id$ | |
2 | // | |
3 | // Emcal jet finder task. | |
4 | // | |
5 | // Authors: C.Loizides, S.Aiola | |
6 | ||
7 | #include <TChain.h> | |
8 | #include <TClonesArray.h> | |
9 | #include <TList.h> | |
10 | #include <TLorentzVector.h> | |
11 | ||
12 | #include "AliAnalysisManager.h" | |
13 | #include "AliCentrality.h" | |
14 | #include "AliEMCALGeometry.h" | |
15 | #include "AliEmcalJet.h" | |
16 | #include "AliFJWrapper.h" | |
17 | #include "AliVCluster.h" | |
18 | #include "AliVParticle.h" | |
19 | #include "AliVEvent.h" | |
20 | #include "AliMCEvent.h" | |
21 | ||
22 | #include "AliEmcalJetTask.h" | |
23 | ||
24 | ClassImp(AliEmcalJetTask) | |
25 | ||
26 | //________________________________________________________________________ | |
27 | AliEmcalJetTask::AliEmcalJetTask() : | |
28 | AliAnalysisTaskSE("AliEmcalJetTask"), | |
29 | fTracksName("Tracks"), | |
30 | fCaloName("CaloClusters"), | |
31 | fJetsName("Jets"), | |
32 | fMC(kFALSE), | |
33 | fAlgo(1), | |
34 | fRadius(0.4), | |
35 | fType(0), | |
36 | fMinJetTrackPt(0.15), | |
37 | fMinJetClusPt(0.15), | |
38 | fMinJetArea(0.01), | |
39 | fMinJetPt(1.0), | |
40 | fJets(0), | |
41 | fEvent(0) | |
42 | { | |
43 | // Default constructor. | |
44 | } | |
45 | ||
46 | //________________________________________________________________________ | |
47 | AliEmcalJetTask::AliEmcalJetTask(const char *name) : | |
48 | AliAnalysisTaskSE(name), | |
49 | fTracksName("Tracks"), | |
50 | fCaloName("CaloClusters"), | |
51 | fJetsName("Jets"), | |
52 | fMC(kFALSE), | |
53 | fAlgo(1), | |
54 | fRadius(0.4), | |
55 | fType(0), | |
56 | fMinJetTrackPt(0.15), | |
57 | fMinJetClusPt(0.15), | |
58 | fMinJetArea(0.01), | |
59 | fMinJetPt(1.0), | |
60 | fJets(0), | |
61 | fEvent(0) | |
62 | { | |
63 | // Standard constructor. | |
64 | ||
65 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex."; | |
66 | } | |
67 | ||
68 | //________________________________________________________________________ | |
69 | AliEmcalJetTask::~AliEmcalJetTask() | |
70 | { | |
71 | // Destructor | |
72 | } | |
73 | ||
74 | //________________________________________________________________________ | |
75 | void AliEmcalJetTask::UserCreateOutputObjects() | |
76 | { | |
77 | // Create user objects. | |
78 | ||
79 | fJets = new TClonesArray("AliEmcalJet"); | |
80 | fJets->SetName(fJetsName); | |
81 | } | |
82 | ||
83 | //________________________________________________________________________ | |
84 | void AliEmcalJetTask::UserExec(Option_t *) | |
85 | { | |
86 | // Main loop, called for each event. | |
87 | ||
88 | // get the event | |
89 | fEvent = InputEvent(); | |
90 | ||
91 | if (!fEvent) { | |
92 | AliError("Could not retrieve event! Returning"); | |
93 | return; | |
94 | } | |
95 | ||
96 | // add jets to event if not yet there | |
97 | if (!(fEvent->FindListObject(fJetsName))) | |
98 | fEvent->AddObject(fJets); | |
99 | ||
100 | // delete jet output | |
101 | fJets->Delete(); | |
102 | ||
103 | if (!fEvent) { | |
104 | AliError(Form("Could not get the event! fMC = %d", fMC)); | |
105 | return; | |
106 | } | |
107 | ||
108 | // get input collections | |
109 | TClonesArray *tracks = 0; | |
110 | TClonesArray *clus = 0; | |
111 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); | |
112 | ||
113 | if ((fType == 0) || (fType == 1)) { | |
114 | if (fTracksName == "Tracks") | |
115 | am->LoadBranch("Tracks"); | |
116 | tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName)); | |
117 | if (!tracks) { | |
118 | AliError(Form("Pointer to tracks %s == 0", fTracksName.Data())); | |
119 | return; | |
120 | } | |
121 | } | |
122 | if ((fType == 0) || (fType == 2)) { | |
123 | if (fCaloName == "CaloClusters") | |
124 | am->LoadBranch("CaloClusters"); | |
125 | clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName)); | |
126 | if (!clus) { | |
127 | AliError(Form("Pointer to clus %s == 0", fCaloName.Data())); | |
128 | return; | |
129 | } | |
130 | } | |
131 | ||
132 | // get centrality | |
133 | Float_t cent = -1; | |
134 | AliCentrality *centrality = fEvent->GetCentrality(); | |
135 | if (centrality) | |
136 | cent = centrality->GetCentralityPercentile("V0M"); | |
137 | else | |
138 | cent = 99; // probably pp data | |
139 | if (cent < 0) { | |
140 | AliWarning(Form("Centrality negative: %f, assuming 99", cent)); | |
141 | cent = 99; | |
142 | } | |
143 | ||
144 | FindJets(tracks, clus, fAlgo, fRadius, cent); | |
145 | } | |
146 | ||
147 | //________________________________________________________________________ | |
148 | void AliEmcalJetTask::Terminate(Option_t *) | |
149 | { | |
150 | // Called once at the end of the analysis. | |
151 | } | |
152 | ||
153 | //________________________________________________________________________ | |
154 | void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/) | |
155 | { | |
156 | // Find jets. | |
157 | ||
158 | if (!tracks && !clus) | |
159 | return; | |
160 | ||
161 | TString name("kt"); | |
162 | fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm); | |
163 | if (algo>=1) { | |
164 | name = "antikt"; | |
165 | jalgo = fastjet::antikt_algorithm; | |
166 | } | |
167 | ||
168 | AliFJWrapper fjw(name, name); | |
169 | fjw.SetAreaType(fastjet::active_area_explicit_ghosts); | |
170 | fjw.SetR(radius); | |
171 | fjw.SetAlgorithm(jalgo); | |
172 | fjw.SetMaxRap(0.9); | |
173 | fjw.Clear(); | |
174 | ||
175 | if (tracks) { | |
176 | const Int_t Ntracks = tracks->GetEntries(); | |
177 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { | |
178 | AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks)); | |
179 | if (!t) | |
180 | continue; | |
181 | ||
182 | if (t->Pt() < fMinJetTrackPt) | |
183 | continue; | |
184 | ||
185 | Int_t index = 1; | |
186 | if(fMC && t->Charge() == 0) | |
187 | index =- 1; | |
188 | ||
189 | fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), (iTracks + 123) * index); | |
190 | } | |
191 | } | |
192 | ||
193 | if (clus) { | |
194 | Double_t vertex[3] = {0, 0, 0}; | |
195 | fEvent->GetPrimaryVertex()->GetXYZ(vertex); | |
196 | const Int_t Nclus = clus->GetEntries(); | |
197 | for (Int_t iClus = 0; iClus < Nclus; ++iClus) { | |
198 | AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus)); | |
199 | if (!c) | |
200 | continue; | |
201 | if (!c->IsEMCAL()) | |
202 | continue; | |
203 | TLorentzVector nPart; | |
204 | c->GetMomentum(nPart, vertex); | |
205 | Double_t et = nPart.Pt(); | |
206 | if (et < fMinJetClusPt) | |
207 | continue; | |
208 | fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 123); | |
209 | } | |
210 | } | |
211 | ||
212 | // run jet finder | |
213 | fjw.Run(); | |
214 | ||
215 | // get geometry | |
216 | AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); | |
217 | if (!geom) { | |
218 | AliFatal("Can not create geometry"); | |
219 | return; | |
220 | } | |
221 | ||
222 | std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets(); | |
223 | for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) { | |
224 | if (jets_incl[ij].perp()<fMinJetPt) | |
225 | continue; | |
226 | if (fjw.GetJetArea(ij)<fMinJetArea) | |
227 | continue; | |
228 | AliEmcalJet *jet = new ((*fJets)[jetCount]) | |
229 | AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m()); | |
230 | Double_t vertex[3] = {0, 0, 0}; | |
231 | fEvent->GetPrimaryVertex()->GetXYZ(vertex); | |
232 | vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij); | |
233 | Int_t nt = 0; | |
234 | Int_t nc = 0; | |
235 | Double_t neutralE = 0; | |
236 | Double_t maxTrack = 0; | |
237 | Double_t maxCluster = 0; | |
238 | Int_t gall = 0; | |
239 | Int_t gemc = 0; | |
240 | ||
241 | jet->SetNumberOfTracks(constituents.size()); | |
242 | jet->SetNumberOfClusters(constituents.size()); | |
243 | ||
244 | for(UInt_t ic = 0; ic < constituents.size(); ++ic) { | |
245 | Int_t uid = constituents[ic].user_index(); | |
246 | ||
247 | if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle | |
248 | ++gall; | |
249 | Double_t gphi = constituents[ic].phi() * TMath::RadToDeg(); | |
250 | Double_t geta = constituents[ic].eta(); | |
251 | if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) && | |
252 | (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax())) | |
253 | ++gemc; | |
254 | continue; | |
255 | } else if (fMC || uid >= 0) { | |
256 | Int_t tid = TMath::Abs(uid) - 123; | |
257 | AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid)); | |
258 | if (t) { | |
259 | if (uid >= 0) | |
260 | neutralE += t->P(); | |
261 | if (t->Pt() > maxTrack) | |
262 | maxTrack = t->Pt(); | |
263 | jet->AddTrackAt(tid, nt); | |
264 | nt++; | |
265 | } | |
266 | } else { | |
267 | Int_t cid = TMath::Abs(uid) - 123; | |
268 | AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid)); | |
269 | if (c) { | |
270 | TLorentzVector nP; | |
271 | c->GetMomentum(nP, vertex); | |
272 | neutralE += nP.P(); | |
273 | if (nP.Pt() > maxCluster) | |
274 | maxCluster = nP.Pt(); | |
275 | jet->AddClusterAt(cid, nc); | |
276 | nc++; | |
277 | } | |
278 | } | |
279 | } | |
280 | jet->SetNumberOfTracks(nt); | |
281 | jet->SetNumberOfClusters(nc); | |
282 | jet->SortConstituents(); | |
283 | jet->SetMaxTrackPt(maxTrack); | |
284 | jet->SetMaxClusterPt(maxCluster); | |
285 | jet->SetNEF(neutralE / jet->E()); | |
286 | jet->SetArea(fjw.GetJetArea(ij)); | |
287 | if (gall > 0) | |
288 | jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall); | |
289 | else | |
290 | jet->SetAreaEmc(-1); | |
291 | if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && | |
292 | (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) && | |
293 | (jet->Eta() > geom->GetArm1EtaMin()) && | |
294 | (jet->Eta() < geom->GetArm1EtaMax())) | |
295 | jet->SetAxisInEmcal(kTRUE); | |
296 | jetCount++; | |
297 | } | |
298 | } |