update from salvatore
[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
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 148void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 149{
150 // Called once at the end of the analysis.
48d874ff 151}
152
153//________________________________________________________________________
914d486c 154void 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}