add task macro
[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
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 24ClassImp(AliEmcalJetTask)
48d874ff 25
48d874ff 26//________________________________________________________________________
914d486c 27AliEmcalJetTask::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
46//________________________________________________________________________
914d486c 47AliEmcalJetTask::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 69AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 70{
71 // Destructor
72}
73
74//________________________________________________________________________
914d486c 75void 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 84void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 85{
86 // Main loop, called for each event.
87
b65fac7a 88 // get the event
89 if (fMC)
90 fEvent = MCEvent();
91 else
92 fEvent = InputEvent();
93
63206072 94 // add jets to event if not yet there
b65fac7a 95 if (!(fEvent->FindListObject(fJetsName)))
96 fEvent->AddObject(fJets);
48d874ff 97
63206072 98 // delete jet output
99 fJets->Delete();
100
b65fac7a 101 if (!fEvent) {
102 AliError(Form("Could not get the event! fMC = %d", fMC));
103 return;
104 }
105
63206072 106 // get input collections
48d874ff 107 TClonesArray *tracks = 0;
108 TClonesArray *clus = 0;
63206072 109 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0 110
b65fac7a 111 if ((fType == 0) || (fType == 1)) {
7efbea04 112 if (fTracksName == "Tracks")
48d874ff 113 am->LoadBranch("Tracks");
b65fac7a 114 tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
48d874ff 115 if (!tracks) {
653159b9 116 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
48d874ff 117 return;
118 }
119 }
b65fac7a 120 if ((fType == 0) || (fType == 2)) {
48d874ff 121 if (fCaloName == "CaloClusters")
122 am->LoadBranch("CaloClusters");
b65fac7a 123 clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
48d874ff 124 if (!clus) {
653159b9 125 AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
48d874ff 126 return;
127 }
48d874ff 128 }
63206072 129
130 // get centrality
131 Float_t cent = -1;
b65fac7a 132 AliCentrality *centrality = fEvent->GetCentrality();
63206072 133 if (centrality)
134 cent = centrality->GetCentralityPercentile("V0M");
135 else
136 cent=99; // probably pp data
b65fac7a 137 if (cent < 0) {
63206072 138 AliError(Form("Centrality negative: %f", cent));
139 return;
140 }
653159b9 141
6ea168d0 142 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 143}
144
145//________________________________________________________________________
914d486c 146void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 147{
148 // Called once at the end of the analysis.
48d874ff 149}
150
151//________________________________________________________________________
914d486c 152void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 153{
154 // Find jets.
155
d03084cd 156 if (!tracks && !clus)
157 return;
158
48d874ff 159 TString name("kt");
160 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
161 if (algo>=1) {
162 name = "antikt";
163 jalgo = fastjet::antikt_algorithm;
164 }
165
166 AliFJWrapper fjw(name, name);
f5f3c8e9 167 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
48d874ff 168 fjw.SetR(radius);
6ea168d0 169 fjw.SetAlgorithm(jalgo);
48d874ff 170 fjw.SetMaxRap(0.9);
171 fjw.Clear();
172
173 if (tracks) {
174 const Int_t Ntracks = tracks->GetEntries();
175 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
b65fac7a 176 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
48d874ff 177 if (!t)
178 continue;
914d486c 179 if (t->Pt()<fMinJetTrackPt)
180 continue;
b65fac7a 181
182 Int_t index = 1;
183 if(fMC && t->Charge() == 0)
184 index =- 1;
185
186 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), (iTracks + 123) * index);
48d874ff 187 }
188 }
6ea168d0 189
48d874ff 190 if (clus) {
191 Double_t vertex[3] = {0, 0, 0};
b65fac7a 192 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 193 const Int_t Nclus = clus->GetEntries();
914d486c 194 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 195 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
d03084cd 196 if (!c)
197 continue;
48d874ff 198 if (!c->IsEMCAL())
199 continue;
200 TLorentzVector nPart;
201 c->GetMomentum(nPart, vertex);
63206072 202 Double_t et = nPart.Pt();
b65fac7a 203 if (et < fMinJetClusPt)
914d486c 204 continue;
b65fac7a 205 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 123);
48d874ff 206 }
207 }
208
7efbea04 209 // run jet finder
48d874ff 210 fjw.Run();
f5f3c8e9 211
212 // get geometry
213 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
214 if (!geom) {
215 AliFatal("Can not create geometry");
216 return;
217 }
6ea168d0 218
48d874ff 219 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a 220 for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9 221 if (jets_incl[ij].perp()<fMinJetPt)
222 continue;
223 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 224 continue;
914d486c 225 AliEmcalJet *jet = new ((*fJets)[jetCount])
226 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0 227 Double_t vertex[3] = {0, 0, 0};
b65fac7a 228 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 229 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
b65fac7a 230 Int_t nt = 0;
231 Int_t nc = 0;
f5f3c8e9 232 Double_t neutralE = 0;
233 Double_t maxTrack = 0;
234 Double_t maxCluster = 0;
235 Int_t gall = 0;
236 Int_t gemc = 0;
b65fac7a 237
f5f3c8e9 238 jet->SetNumberOfTracks(constituents.size());
239 jet->SetNumberOfClusters(constituents.size());
b65fac7a 240
241 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff 242 Int_t uid = constituents[ic].user_index();
b65fac7a 243
244 if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
f5f3c8e9 245 ++gall;
b65fac7a 246 Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
f5f3c8e9 247 Double_t geta = constituents[ic].eta();
b65fac7a 248 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
249 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
f5f3c8e9 250 ++gemc;
251 continue;
b65fac7a 252 } else if (fMC || uid >= 0) {
253 Int_t tid = TMath::Abs(uid) - 123;
254 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
d03084cd 255 if (t) {
b65fac7a 256 if (uid >= 0)
257 neutralE += t->P();
258 if (t->Pt() > maxTrack)
259 maxTrack = t->Pt();
260 jet->AddTrackAt(tid, nt);
d03084cd 261 nt++;
262 }
6ea168d0 263 } else {
b65fac7a 264 Int_t cid = TMath::Abs(uid) - 123;
265 AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
266 if (c) {
267 TLorentzVector nP;
268 c->GetMomentum(nP, vertex);
269 neutralE += nP.P();
270 if (nP.Pt() > maxCluster)
271 maxCluster = nP.Pt();
272 jet->AddClusterAt(cid, nc);
273 nc++;
274 }
48d874ff 275 }
276 }
35789a2d 277 jet->SetNumberOfTracks(nt);
278 jet->SetNumberOfClusters(nc);
f5f3c8e9 279 jet->SortConstituents();
6ea168d0 280 jet->SetMaxTrackPt(maxTrack);
281 jet->SetMaxClusterPt(maxCluster);
b65fac7a 282 jet->SetNEF(neutralE / jet->E());
f5f3c8e9 283 jet->SetArea(fjw.GetJetArea(ij));
b65fac7a 284 if (gall > 0)
285 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 286 else
287 jet->SetAreaEmc(-1);
b65fac7a 288 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
289 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
290 (jet->Eta() > geom->GetArm1EtaMin()) &&
291 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 292 jet->SetAxisInEmcal(kTRUE);
48d874ff 293 jetCount++;
294 }
295}