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