updates for train
[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"
111318e8 20#include "AliESDEvent.h"
b65fac7a 21#include "AliMCEvent.h"
22
23#include "AliEmcalJetTask.h"
48d874ff 24
914d486c 25ClassImp(AliEmcalJetTask)
48d874ff 26
48d874ff 27//________________________________________________________________________
914d486c 28AliEmcalJetTask::AliEmcalJetTask() :
29 AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0 30 fTracksName("Tracks"),
31 fCaloName("CaloClusters"),
32 fJetsName("Jets"),
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"),
48d874ff 52 fAlgo(1),
53 fRadius(0.4),
54 fType(0),
6ea168d0 55 fMinJetTrackPt(0.15),
914d486c 56 fMinJetClusPt(0.15),
f5f3c8e9 57 fMinJetArea(0.01),
58 fMinJetPt(1.0),
b65fac7a 59 fJets(0),
60 fEvent(0)
48d874ff 61{
62 // Standard constructor.
7efbea04 63
7efbea04 64 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff 65}
66
67//________________________________________________________________________
914d486c 68AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 69{
70 // Destructor
71}
72
73//________________________________________________________________________
914d486c 74void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff 75{
76 // Create user objects.
77
914d486c 78 fJets = new TClonesArray("AliEmcalJet");
48d874ff 79 fJets->SetName(fJetsName);
48d874ff 80}
81
111318e8 82//_____________________________________________________
83TString AliEmcalJetTask::GetBeamType()
84{
85 // Get beam type : pp-AA-pA
86 // ESDs have it directly, AODs get it from hardcoded run number ranges
87
88 TString beamType;
89
90 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
91 if (esd) {
92 const AliESDRun *run = esd->GetESDRun();
93 beamType = run->GetBeamType();
94 }
95 else
96 {
97 Int_t runNumber = fEvent->GetRunNumber();
98 if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
99 (runNumber >= 166529 && runNumber <= 170593)) // LHC11h
100 {
101 beamType = "A-A";
102 }
103 else
104 {
105 beamType = "p-p";
106 }
107 }
108
109 return beamType;
110}
111
112
48d874ff 113//________________________________________________________________________
914d486c 114void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 115{
116 // Main loop, called for each event.
117
b65fac7a 118 // get the event
1772fe65 119 fEvent = InputEvent();
120
121 if (!fEvent) {
02c7e943 122 AliError(Form("%s: Could not retrieve event! Returning", GetName()));
1772fe65 123 return;
124 }
b65fac7a 125
63206072 126 // add jets to event if not yet there
b65fac7a 127 if (!(fEvent->FindListObject(fJetsName)))
128 fEvent->AddObject(fJets);
48d874ff 129
63206072 130 // delete jet output
131 fJets->Delete();
132
63206072 133 // get input collections
48d874ff 134 TClonesArray *tracks = 0;
135 TClonesArray *clus = 0;
63206072 136 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0 137
b65fac7a 138 if ((fType == 0) || (fType == 1)) {
7efbea04 139 if (fTracksName == "Tracks")
48d874ff 140 am->LoadBranch("Tracks");
b65fac7a 141 tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
48d874ff 142 if (!tracks) {
02c7e943 143 AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
48d874ff 144 return;
145 }
146 }
b65fac7a 147 if ((fType == 0) || (fType == 2)) {
48d874ff 148 if (fCaloName == "CaloClusters")
149 am->LoadBranch("CaloClusters");
b65fac7a 150 clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
48d874ff 151 if (!clus) {
02c7e943 152 AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
48d874ff 153 return;
154 }
48d874ff 155 }
111318e8 156
157 // get centrality
158 Double_t cent = 99;
159
160 if (GetBeamType() == "A-A") {
161 AliCentrality *centrality = InputEvent()->GetCentrality();
162
163 if (centrality)
164 cent = centrality->GetCentralityPercentile("V0M");
165 else
166 cent = 99; // probably pp data
167
168 if (cent < 0) {
169 AliWarning(Form("Centrality negative: %f, assuming 99", cent));
170 cent = 99;
171 }
63206072 172 }
653159b9 173
6ea168d0 174 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 175}
176
177//________________________________________________________________________
914d486c 178void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 179{
180 // Called once at the end of the analysis.
48d874ff 181}
182
183//________________________________________________________________________
914d486c 184void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 185{
186 // Find jets.
187
d03084cd 188 if (!tracks && !clus)
189 return;
190
48d874ff 191 TString name("kt");
192 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
193 if (algo>=1) {
194 name = "antikt";
195 jalgo = fastjet::antikt_algorithm;
196 }
197
198 AliFJWrapper fjw(name, name);
f5f3c8e9 199 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
48d874ff 200 fjw.SetR(radius);
6ea168d0 201 fjw.SetAlgorithm(jalgo);
48d874ff 202 fjw.SetMaxRap(0.9);
203 fjw.Clear();
204
205 if (tracks) {
206 const Int_t Ntracks = tracks->GetEntries();
207 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1772fe65 208 AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
48d874ff 209 if (!t)
210 continue;
1772fe65 211
212 if (t->Pt() < fMinJetTrackPt)
914d486c 213 continue;
a3c8c8c8 214 // offset of 100 for consistency with cluster ids
215 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
48d874ff 216 }
217 }
6ea168d0 218
48d874ff 219 if (clus) {
220 Double_t vertex[3] = {0, 0, 0};
b65fac7a 221 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 222 const Int_t Nclus = clus->GetEntries();
914d486c 223 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 224 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
d03084cd 225 if (!c)
226 continue;
48d874ff 227 if (!c->IsEMCAL())
228 continue;
229 TLorentzVector nPart;
230 c->GetMomentum(nPart, vertex);
111318e8 231 if (nPart.Pt() < fMinJetClusPt)
914d486c 232 continue;
a3c8c8c8 233 // offset of 100 to skip ghost particles uid = -1
234 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);
48d874ff 235 }
236 }
7efbea04 237 // run jet finder
48d874ff 238 fjw.Run();
f5f3c8e9 239
240 // get geometry
241 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
242 if (!geom) {
02c7e943 243 AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9 244 return;
245 }
6ea168d0 246
48d874ff 247 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a 248 for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9 249 if (jets_incl[ij].perp()<fMinJetPt)
250 continue;
251 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 252 continue;
914d486c 253 AliEmcalJet *jet = new ((*fJets)[jetCount])
254 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0 255 Double_t vertex[3] = {0, 0, 0};
b65fac7a 256 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 257 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
b65fac7a 258 Int_t nt = 0;
259 Int_t nc = 0;
f5f3c8e9 260 Double_t neutralE = 0;
111318e8 261 Double_t maxCh = 0;
262 Double_t maxNe = 0;
f5f3c8e9 263 Int_t gall = 0;
264 Int_t gemc = 0;
111318e8 265 Int_t ncharged = 0;
266 Int_t nneutral = 0;
02c7e943 267 Double_t mcpt = 0;
b65fac7a 268
f5f3c8e9 269 jet->SetNumberOfTracks(constituents.size());
270 jet->SetNumberOfClusters(constituents.size());
b65fac7a 271 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff 272 Int_t uid = constituents[ic].user_index();
b65fac7a 273
274 if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
f5f3c8e9 275 ++gall;
b65fac7a 276 Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
f5f3c8e9 277 Double_t geta = constituents[ic].eta();
b65fac7a 278 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
279 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
f5f3c8e9 280 ++gemc;
281 continue;
02c7e943 282 } else if ((uid > 0) && tracks) { // track constituent
111318e8 283 Int_t tid = uid - 100;
b65fac7a 284 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
d03084cd 285 if (t) {
111318e8 286 if (t->Charge() == 0) {
b65fac7a 287 neutralE += t->P();
111318e8 288 ++nneutral;
289 if (t->Pt() > maxNe)
290 maxNe = t->Pt();
291 } else {
292 ++ncharged;
293 if (t->Pt() > maxCh)
294 maxCh = t->Pt();
295 }
296
297 if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
02c7e943 298 mcpt += t->Pt();
111318e8 299
b65fac7a 300 jet->AddTrackAt(tid, nt);
111318e8 301 ++nt;
d03084cd 302 }
02c7e943 303 } else if (clus) { // cluster constituent
111318e8 304 Int_t cid = -uid - 100;
b65fac7a 305 AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
306 if (c) {
307 TLorentzVector nP;
308 c->GetMomentum(nP, vertex);
309 neutralE += nP.P();
111318e8 310 if (nP.Pt() > maxNe)
311 maxNe = nP.Pt();
312
313 if (c->Chi2() == 100) // MC particle
02c7e943 314 mcpt += nP.Pt();
111318e8 315
b65fac7a 316 jet->AddClusterAt(cid, nc);
111318e8 317 ++nc;
318 ++nneutral;
b65fac7a 319 }
02c7e943 320 } else {
321 AliError(Form("%s: No logical way to end up here.", GetName()));
322 continue;
48d874ff 323 }
324 }
35789a2d 325 jet->SetNumberOfTracks(nt);
326 jet->SetNumberOfClusters(nc);
f5f3c8e9 327 jet->SortConstituents();
111318e8 328 jet->SetMaxNeutralPt(maxNe);
329 jet->SetMaxChargedPt(maxCh);
b65fac7a 330 jet->SetNEF(neutralE / jet->E());
f5f3c8e9 331 jet->SetArea(fjw.GetJetArea(ij));
111318e8 332 jet->SetNumberOfCharged(ncharged);
333 jet->SetNumberOfNeutrals(nneutral);
02c7e943 334 jet->SetMCPt(mcpt);
b65fac7a 335 if (gall > 0)
336 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 337 else
338 jet->SetAreaEmc(-1);
b65fac7a 339 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65 340 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
341 (jet->Eta() > geom->GetArm1EtaMin()) &&
342 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 343 jet->SetAxisInEmcal(kTRUE);
48d874ff 344 jetCount++;
345 }
346}