vector area and compare fct.
[u/mrichter/AliRoot.git] / PWGGA / 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
63206072 113 // delete jet output
114 fJets->Delete();
115
04905231 116 FindJets();
48d874ff 117}
118
119//________________________________________________________________________
914d486c 120void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 121{
122 // Called once at the end of the analysis.
48d874ff 123}
124
125//________________________________________________________________________
04905231 126void 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//________________________________________________________________________
396Bool_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}