Changes by Salvatore
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx
CommitLineData
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 27ClassImp(AliEmcalJetTask)
48d874ff 28
48d874ff 29//________________________________________________________________________
914d486c 30AliEmcalJetTask::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 58AliEmcalJetTask::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 88AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 89{
90 // Destructor
91}
92
93//________________________________________________________________________
914d486c 94void 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 103void 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
04905231 113 FindJets();
48d874ff 114}
115
116//________________________________________________________________________
914d486c 117void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 118{
119 // Called once at the end of the analysis.
48d874ff 120}
121
122//________________________________________________________________________
04905231 123void AliEmcalJetTask::FindJets()
48d874ff 124{
125 // Find jets.
126
04905231 127 if (!fTracks && !fClus)
d03084cd 128 return;
129
48d874ff 130 TString name("kt");
131 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
04905231 132 if (fAlgo>=1) {
48d874ff 133 name = "antikt";
134 jalgo = fastjet::antikt_algorithm;
135 }
136
04905231 137 // setup fj wrapper
48d874ff 138 AliFJWrapper fjw(name, name);
f5f3c8e9 139 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
1f9c287f 140 fjw.SetGhostArea(fGhostArea);
04905231 141 fjw.SetR(fRadius);
6ea168d0 142 fjw.SetAlgorithm(jalgo);
04905231 143 fjw.SetMaxRap(1);
48d874ff 144 fjw.Clear();
145
04905231 146 // get primary vertex
147 Double_t vertex[3] = {0, 0, 0};
148 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
149
da67b88a 150 if ((fIsMcPart || fType == 0 || fType == 1) && fTracks) {
04905231 151 const Int_t Ntracks = fTracks->GetEntries();
48d874ff 152 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
04905231 153 AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
48d874ff 154 if (!t)
155 continue;
5d950148 156 if (fIsMcPart && fType == 1 && t->Charge() == 0)
da67b88a 157 continue;
5d950148 158 if (fIsMcPart && fType == 2 && t->Charge() != 0)
da67b88a 159 continue;
1772fe65 160 if (t->Pt() < fMinJetTrackPt)
914d486c 161 continue;
04905231 162 Double_t eta = t->Eta();
163 Double_t phi = t->Phi();
ddf74a88 164 if ((eta<fEtaMin) || (eta>fEtaMax) ||
165 (phi<fPhiMin) || (phi>fPhiMax))
04905231 166 continue;
da67b88a 167
a3c8c8c8 168 // offset of 100 for consistency with cluster ids
169 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
48d874ff 170 }
171 }
6ea168d0 172
da67b88a 173 if ((fType == 0 || fType == 2) && fClus) {
174 const Int_t Nclus = fClus->GetEntries();
175 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
176 AliVCluster *c = 0;
177 Double_t cEta=0,cPhi=0,cPt=0;
178 Double_t cPx=0,cPy=0,cPz=0;
179 if (fIsEmcPart) {
180 AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
181 if (!ep)
182 continue;
183 c = ep->GetCluster();
184 if (!c)
04905231 185 continue;
da67b88a 186 cEta = ep->Eta();
187 cPhi = ep->Phi();
188 cPt = ep->Pt();
189 cPx = ep->Px();
190 cPy = ep->Py();
191 cPz = ep->Pz();
192 } else {
193 c = static_cast<AliVCluster*>(fClus->At(iClus));
194 if (!c)
195 continue;
196 TLorentzVector nP;
197 c->GetMomentum(nP, vertex);
198 cEta = nP.Eta();
199 cPhi = nP.Phi();
200 cPt = nP.Pt();
201 cPx = nP.Px();
202 cPy = nP.Py();
203 cPz = nP.Pz();
04905231 204 }
da67b88a 205 if (!c->IsEMCAL())
206 continue;
207 if (cPt < fMinJetClusPt)
208 continue;
ddf74a88 209 if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
210 (cPhi<fPhiMin) || (cPhi>fPhiMax))
da67b88a 211 continue;
212 // offset of 100 to skip ghost particles uid = -1
213 fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
48d874ff 214 }
215 }
da67b88a 216
7efbea04 217 // run jet finder
48d874ff 218 fjw.Run();
f5f3c8e9 219
220 // get geometry
221 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
222 if (!geom) {
02c7e943 223 AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9 224 return;
225 }
04905231 226
227 // loop over fastjet jets
48d874ff 228 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a 229 for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9 230 if (jets_incl[ij].perp()<fMinJetPt)
231 continue;
232 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 233 continue;
04905231 234
914d486c 235 AliEmcalJet *jet = new ((*fJets)[jetCount])
236 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
04905231 237
238 // loop over constituents
c64cb1f6 239 std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231 240 jet->SetNumberOfTracks(constituents.size());
241 jet->SetNumberOfClusters(constituents.size());
242
b65fac7a 243 Int_t nt = 0;
244 Int_t nc = 0;
f5f3c8e9 245 Double_t neutralE = 0;
111318e8 246 Double_t maxCh = 0;
247 Double_t maxNe = 0;
f5f3c8e9 248 Int_t gall = 0;
249 Int_t gemc = 0;
04905231 250 Int_t cemc = 0;
111318e8 251 Int_t ncharged = 0;
252 Int_t nneutral = 0;
02c7e943 253 Double_t mcpt = 0;
04905231 254 Double_t emcpt = 0;
b65fac7a 255
b65fac7a 256 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff 257 Int_t uid = constituents[ic].user_index();
b65fac7a 258
04905231 259 if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
f5f3c8e9 260 ++gall;
04905231 261 Double_t gphi = constituents[ic].phi();
262 if (gphi<0)
263 gphi += TMath::TwoPi();
264 gphi *= TMath::RadToDeg();
f5f3c8e9 265 Double_t geta = constituents[ic].eta();
b65fac7a 266 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
267 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231 268 ++gemc;
269 } else if ((uid > 0) && fTracks) { // track constituent
111318e8 270 Int_t tid = uid - 100;
04905231 271 AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
272 if (!t)
273 continue;
274 Double_t cEta = t->Eta();
275 Double_t cPhi = t->Phi();
276 Double_t cPt = t->Pt();
277 Double_t cP = t->P();
278 if (t->Charge() == 0) {
279 neutralE += cP;
280 ++nneutral;
281 if (cPt > maxNe)
282 maxNe = cPt;
283 } else {
284 ++ncharged;
285 if (cPt > maxCh)
286 maxCh = cPt;
287 }
288
289 if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
290 mcpt += cPt;
291
292 if (cPhi<0)
293 cPhi += TMath::TwoPi();
294 cPhi *= TMath::RadToDeg();
295 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
296 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
297 emcpt += cPt;
298 ++cemc;
d03084cd 299 }
04905231 300
301 jet->AddTrackAt(tid, nt);
302 ++nt;
303 } else if (fClus) { // cluster constituent
111318e8 304 Int_t cid = -uid - 100;
04905231 305 AliVCluster *c = 0;
306 Double_t cEta=0,cPhi=0,cPt=0,cP=0;
307 if (fIsEmcPart) {
308 AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
309 if (!ep)
310 continue;
311 c = ep->GetCluster();
312 if (!c)
313 continue;
314 cEta = ep->Eta();
315 cPhi = ep->Phi();
316 cPt = ep->Pt();
317 cP = ep->P();
318 } else {
319 c = static_cast<AliVCluster*>(fClus->At(cid));
320 if (!c)
321 continue;
322 TLorentzVector nP;
323 c->GetMomentum(nP, vertex);
324 cEta = nP.Eta();
325 cPhi = nP.Phi();
326 cPt = nP.Pt();
327 cP = nP.P();
328 }
329
330 neutralE += cP;
331 if (cPt > maxNe)
332 maxNe = cPt;
333
334 if (c->Chi2() == 100) // MC particle
335 mcpt += cPt;
336
337 if (cPhi<0)
338 cPhi += TMath::TwoPi();
339 cPhi *= TMath::RadToDeg();
340 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
341 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
342 emcpt += cPt;
343 ++cemc;
344 }
345
346 jet->AddClusterAt(cid, nc);
347 ++nc;
348 ++nneutral;
02c7e943 349 } else {
350 AliError(Form("%s: No logical way to end up here.", GetName()));
351 continue;
48d874ff 352 }
04905231 353 } /* end of constituent loop */
354
35789a2d 355 jet->SetNumberOfTracks(nt);
356 jet->SetNumberOfClusters(nc);
f5f3c8e9 357 jet->SortConstituents();
111318e8 358 jet->SetMaxNeutralPt(maxNe);
359 jet->SetMaxChargedPt(maxCh);
b65fac7a 360 jet->SetNEF(neutralE / jet->E());
a88b5c99 361 fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
362 jet->SetArea(area.perp());
363 jet->SetAreaEta(area.eta());
364 jet->SetAreaPhi(area.phi());
111318e8 365 jet->SetNumberOfCharged(ncharged);
366 jet->SetNumberOfNeutrals(nneutral);
02c7e943 367 jet->SetMCPt(mcpt);
04905231 368 jet->SetNEmc(cemc);
369 jet->SetPtEmc(emcpt);
370
b65fac7a 371 if (gall > 0)
372 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 373 else
374 jet->SetAreaEmc(-1);
b65fac7a 375 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65 376 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
377 (jet->Eta() > geom->GetArm1EtaMin()) &&
378 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 379 jet->SetAxisInEmcal(kTRUE);
48d874ff 380 jetCount++;
381 }
a88b5c99 382 fJets->Sort();
48d874ff 383}
04905231 384
385//________________________________________________________________________
386Bool_t AliEmcalJetTask::DoInit()
387{
388 // Init. Return true if successful.
389
390 // get input collections
391 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
392
393 // get the event
394 fEvent = InputEvent();
395 if (!fEvent) {
396 AliError(Form("%s: Could not retrieve event! Returning", GetName()));
397 return 0;
398 }
399
400 // add jets to event if not yet there
71c1dd76 401 fJets->Delete();
04905231 402 if (!(fEvent->FindListObject(fJetsName)))
403 fEvent->AddObject(fJets);
404 else {
f4b49fa6 405 AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
04905231 406 return 0;
407 }
408
da67b88a 409 if (fTracksName == "Tracks")
410 am->LoadBranch("Tracks");
411 if (!fTracks && !fTracksName.IsNull()) {
412 fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
413 if (!fTracks) {
414 AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
415 return 0;
04905231 416 }
da67b88a 417 }
418 if (fTracks) {
419 TClass cls(fTracks->GetClass()->GetName());
04905231 420 if (cls.InheritsFrom("AliMCParticle"))
421 fIsMcPart = 1;
422 }
da67b88a 423
424 if (fCaloName == "CaloClusters")
425 am->LoadBranch("CaloClusters");
426 if (!fClus && !fCaloName.IsNull()) {
427 fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
428 if (!fClus) {
429 AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
430 return 0;
04905231 431 }
da67b88a 432 }
433 if (fClus) {
434 TClass cls(fClus->GetClass()->GetName());
04905231 435 if (cls.InheritsFrom("AliEmcalParticle"))
436 fIsEmcPart = 1;
437 }
438
439 return 1;
440}