New histograms added for Raws:
[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) {
122 AliError("Could not retrieve event! Returning");
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) {
653159b9 143 AliError(Form("Pointer to tracks %s == 0", 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) {
653159b9 152 AliError(Form("Pointer to clus %s == 0", 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 }
237
7efbea04 238 // run jet finder
48d874ff 239 fjw.Run();
f5f3c8e9 240
241 // get geometry
242 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
243 if (!geom) {
244 AliFatal("Can not create geometry");
245 return;
246 }
6ea168d0 247
48d874ff 248 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
b65fac7a 249 for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9 250 if (jets_incl[ij].perp()<fMinJetPt)
251 continue;
252 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 253 continue;
914d486c 254 AliEmcalJet *jet = new ((*fJets)[jetCount])
255 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0 256 Double_t vertex[3] = {0, 0, 0};
b65fac7a 257 fEvent->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 258 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
b65fac7a 259 Int_t nt = 0;
260 Int_t nc = 0;
f5f3c8e9 261 Double_t neutralE = 0;
111318e8 262 Double_t maxCh = 0;
263 Double_t maxNe = 0;
f5f3c8e9 264 Int_t gall = 0;
265 Int_t gemc = 0;
111318e8 266 Int_t ncharged = 0;
267 Int_t nneutral = 0;
268 Double_t MCpt = 0;
b65fac7a 269
f5f3c8e9 270 jet->SetNumberOfTracks(constituents.size());
271 jet->SetNumberOfClusters(constituents.size());
b65fac7a 272 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff 273 Int_t uid = constituents[ic].user_index();
b65fac7a 274
275 if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
f5f3c8e9 276 ++gall;
b65fac7a 277 Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
f5f3c8e9 278 Double_t geta = constituents[ic].eta();
b65fac7a 279 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
280 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
f5f3c8e9 281 ++gemc;
282 continue;
a3c8c8c8 283 } else if (uid > 0) { // track constituent
111318e8 284 Int_t tid = uid - 100;
b65fac7a 285 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
d03084cd 286 if (t) {
111318e8 287 if (t->Charge() == 0) {
b65fac7a 288 neutralE += t->P();
111318e8 289 ++nneutral;
290 if (t->Pt() > maxNe)
291 maxNe = t->Pt();
292 } else {
293 ++ncharged;
294 if (t->Pt() > maxCh)
295 maxCh = t->Pt();
296 }
297
298 if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
299 MCpt += t->Pt();
300
b65fac7a 301 jet->AddTrackAt(tid, nt);
111318e8 302 ++nt;
d03084cd 303 }
a3c8c8c8 304 } else { // cluster constituent
111318e8 305 Int_t cid = -uid - 100;
b65fac7a 306 AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
307 if (c) {
308 TLorentzVector nP;
309 c->GetMomentum(nP, vertex);
310 neutralE += nP.P();
111318e8 311 if (nP.Pt() > maxNe)
312 maxNe = nP.Pt();
313
314 if (c->Chi2() == 100) // MC particle
315 MCpt += nP.Pt();
316
b65fac7a 317 jet->AddClusterAt(cid, nc);
111318e8 318 ++nc;
319 ++nneutral;
b65fac7a 320 }
48d874ff 321 }
322 }
35789a2d 323 jet->SetNumberOfTracks(nt);
324 jet->SetNumberOfClusters(nc);
f5f3c8e9 325 jet->SortConstituents();
111318e8 326 jet->SetMaxNeutralPt(maxNe);
327 jet->SetMaxChargedPt(maxCh);
b65fac7a 328 jet->SetNEF(neutralE / jet->E());
f5f3c8e9 329 jet->SetArea(fjw.GetJetArea(ij));
111318e8 330 jet->SetNumberOfCharged(ncharged);
331 jet->SetNumberOfNeutrals(nneutral);
332 jet->SetMCPt(MCpt);
b65fac7a 333 if (gall > 0)
334 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 335 else
336 jet->SetAreaEmc(-1);
b65fac7a 337 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65 338 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
339 (jet->Eta() > geom->GetArm1EtaMin()) &&
340 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 341 jet->SetAxisInEmcal(kTRUE);
48d874ff 342 jetCount++;
343 }
344}