3 // Emcal jet finder task.
5 // Authors: C.Loizides, S.Aiola
8 #include <TClonesArray.h>
10 #include <TLorentzVector.h>
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 "AliESDEvent.h"
21 #include "AliMCEvent.h"
23 #include "AliEmcalJetTask.h"
25 ClassImp(AliEmcalJetTask)
27 //________________________________________________________________________
28 AliEmcalJetTask::AliEmcalJetTask() :
29 AliAnalysisTaskSE("AliEmcalJetTask"),
30 fTracksName("Tracks"),
31 fCaloName("CaloClusters"),
43 // Default constructor.
46 //________________________________________________________________________
47 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
48 AliAnalysisTaskSE(name),
49 fTracksName("Tracks"),
50 fCaloName("CaloClusters"),
62 // Standard constructor.
64 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
67 //________________________________________________________________________
68 AliEmcalJetTask::~AliEmcalJetTask()
73 //________________________________________________________________________
74 void AliEmcalJetTask::UserCreateOutputObjects()
76 // Create user objects.
78 fJets = new TClonesArray("AliEmcalJet");
79 fJets->SetName(fJetsName);
82 //_____________________________________________________
83 TString AliEmcalJetTask::GetBeamType()
85 // Get beam type : pp-AA-pA
86 // ESDs have it directly, AODs get it from hardcoded run number ranges
90 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
92 const AliESDRun *run = esd->GetESDRun();
93 beamType = run->GetBeamType();
97 Int_t runNumber = fEvent->GetRunNumber();
98 if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
99 (runNumber >= 166529 && runNumber <= 170593)) // LHC11h
113 //________________________________________________________________________
114 void AliEmcalJetTask::UserExec(Option_t *)
116 // Main loop, called for each event.
119 fEvent = InputEvent();
122 AliError(Form("%s: Could not retrieve event! Returning", GetName()));
126 // add jets to event if not yet there
127 if (!(fEvent->FindListObject(fJetsName)))
128 fEvent->AddObject(fJets);
133 // get input collections
134 TClonesArray *tracks = 0;
135 TClonesArray *clus = 0;
136 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
138 if ((fType == 0) || (fType == 1)) {
139 if (fTracksName == "Tracks")
140 am->LoadBranch("Tracks");
141 tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
143 AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
147 if ((fType == 0) || (fType == 2)) {
148 if (fCaloName == "CaloClusters")
149 am->LoadBranch("CaloClusters");
150 clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
152 AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
160 if (GetBeamType() == "A-A") {
161 AliCentrality *centrality = InputEvent()->GetCentrality();
164 cent = centrality->GetCentralityPercentile("V0M");
166 cent = 99; // probably pp data
169 AliWarning(Form("Centrality negative: %f, assuming 99", cent));
174 FindJets(tracks, clus, fAlgo, fRadius, cent);
177 //________________________________________________________________________
178 void AliEmcalJetTask::Terminate(Option_t *)
180 // Called once at the end of the analysis.
183 //________________________________________________________________________
184 void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
188 if (!tracks && !clus)
192 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
195 jalgo = fastjet::antikt_algorithm;
198 AliFJWrapper fjw(name, name);
199 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
201 fjw.SetAlgorithm(jalgo);
206 const Int_t Ntracks = tracks->GetEntries();
207 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
208 AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
212 if (t->Pt() < fMinJetTrackPt)
214 // offset of 100 for consistency with cluster ids
215 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
220 Double_t vertex[3] = {0, 0, 0};
221 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
222 const Int_t Nclus = clus->GetEntries();
223 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
224 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
229 TLorentzVector nPart;
230 c->GetMomentum(nPart, vertex);
231 if (nPart.Pt() < fMinJetClusPt)
233 // offset of 100 to skip ghost particles uid = -1
234 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);
242 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
244 AliFatal(Form("%s: Can not create geometry", GetName()));
248 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
249 for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
250 if (jets_incl[ij].perp()<fMinJetPt)
252 if (fjw.GetJetArea(ij)<fMinJetArea)
254 AliEmcalJet *jet = new ((*fJets)[jetCount])
255 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
256 Double_t vertex[3] = {0, 0, 0};
257 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
258 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
261 Double_t neutralE = 0;
270 jet->SetNumberOfTracks(constituents.size());
271 jet->SetNumberOfClusters(constituents.size());
272 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
273 Int_t uid = constituents[ic].user_index();
275 if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
277 Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
278 Double_t geta = constituents[ic].eta();
279 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
280 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
283 } else if ((uid > 0) && tracks) { // track constituent
284 Int_t tid = uid - 100;
285 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
287 if (t->Charge() == 0) {
298 if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
301 jet->AddTrackAt(tid, nt);
304 } else if (clus) { // cluster constituent
305 Int_t cid = -uid - 100;
306 AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
309 c->GetMomentum(nP, vertex);
314 if (c->Chi2() == 100) // MC particle
317 jet->AddClusterAt(cid, nc);
322 AliError(Form("%s: No logical way to end up here.", GetName()));
326 jet->SetNumberOfTracks(nt);
327 jet->SetNumberOfClusters(nc);
328 jet->SortConstituents();
329 jet->SetMaxNeutralPt(maxNe);
330 jet->SetMaxChargedPt(maxCh);
331 jet->SetNEF(neutralE / jet->E());
332 jet->SetArea(fjw.GetJetArea(ij));
333 jet->SetNumberOfCharged(ncharged);
334 jet->SetNumberOfNeutrals(nneutral);
337 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
340 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
341 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
342 (jet->Eta() > geom->GetArm1EtaMin()) &&
343 (jet->Eta() < geom->GetArm1EtaMax()))
344 jet->SetAxisInEmcal(kTRUE);