updates from Salvatore
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliEmcalJetTask.cxx
CommitLineData
48d874ff 1// $Id$
542ac102 2//
3// Emcal jet finder task.
4//
68263112 5// Author: C.Loizides
48d874ff 6
914d486c 7#include "AliEmcalJetTask.h"
6ea168d0 8#include <TChain.h>
48d874ff 9#include <TClonesArray.h>
6ea168d0 10#include <TH1F.h>
11#include <TH2F.h>
12#include <TList.h>
13#include <TLorentzVector.h>
48d874ff 14#include <TParticle.h>
48d874ff 15#include "AliAnalysisManager.h"
6ea168d0 16#include "AliCentrality.h"
f5f3c8e9 17#include "AliEMCALGeometry.h"
18#include "AliEmcalJet.h"
63206072 19#include "AliFJWrapper.h"
35789a2d 20#include "AliVCluster.h"
21#include "AliVTrack.h"
48d874ff 22
914d486c 23ClassImp(AliEmcalJetTask)
48d874ff 24
48d874ff 25//________________________________________________________________________
914d486c 26AliEmcalJetTask::AliEmcalJetTask() :
27 AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0 28 fTracksName("Tracks"),
29 fCaloName("CaloClusters"),
30 fJetsName("Jets"),
31 fAlgo(1),
32 fRadius(0.4),
33 fType(0),
914d486c 34 fMinJetTrackPt(0.15),
35 fMinJetClusPt(0.15),
f5f3c8e9 36 fMinJetArea(0.01),
37 fMinJetPt(1.0),
6ea168d0 38 fJets(0)
39{
40 // Default constructor.
6ea168d0 41}
42
43//________________________________________________________________________
914d486c 44AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072 45 AliAnalysisTaskSE(name),
7efbea04 46 fTracksName("Tracks"),
48d874ff 47 fCaloName("CaloClusters"),
7efbea04 48 fJetsName("Jets"),
48d874ff 49 fAlgo(1),
50 fRadius(0.4),
51 fType(0),
6ea168d0 52 fMinJetTrackPt(0.15),
914d486c 53 fMinJetClusPt(0.15),
f5f3c8e9 54 fMinJetArea(0.01),
55 fMinJetPt(1.0),
7efbea04 56 fJets(0)
48d874ff 57{
58 // Standard constructor.
7efbea04 59
7efbea04 60 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff 61}
62
63//________________________________________________________________________
914d486c 64AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 65{
66 // Destructor
67}
68
69//________________________________________________________________________
914d486c 70void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff 71{
72 // Create user objects.
73
914d486c 74 fJets = new TClonesArray("AliEmcalJet");
48d874ff 75 fJets->SetName(fJetsName);
48d874ff 76}
77
78//________________________________________________________________________
914d486c 79void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 80{
81 // Main loop, called for each event.
82
63206072 83 // add jets to event if not yet there
48d874ff 84 if (!(InputEvent()->FindListObject(fJetsName)))
85 InputEvent()->AddObject(fJets);
86
63206072 87 // delete jet output
88 fJets->Delete();
89
63206072 90 // get input collections
48d874ff 91 TClonesArray *tracks = 0;
92 TClonesArray *clus = 0;
93 TList *l = InputEvent()->GetList();
63206072 94 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0 95
48d874ff 96 if ((fType==0)||(fType==1)) {
7efbea04 97 if (fTracksName == "Tracks")
48d874ff 98 am->LoadBranch("Tracks");
7efbea04 99 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
48d874ff 100 if (!tracks) {
653159b9 101 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
48d874ff 102 return;
103 }
104 }
105 if ((fType==0)||(fType==2)) {
106 if (fCaloName == "CaloClusters")
107 am->LoadBranch("CaloClusters");
108 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
109 if (!clus) {
653159b9 110 AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
48d874ff 111 return;
112 }
48d874ff 113 }
63206072 114
115 // get centrality
116 Float_t cent = -1;
653159b9 117 AliCentrality *centrality = InputEvent()->GetCentrality();
63206072 118 if (centrality)
119 cent = centrality->GetCentralityPercentile("V0M");
120 else
121 cent=99; // probably pp data
122 if (cent<0) {
123 AliError(Form("Centrality negative: %f", cent));
124 return;
125 }
653159b9 126
6ea168d0 127 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 128}
129
130//________________________________________________________________________
914d486c 131void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 132{
133 // Called once at the end of the analysis.
48d874ff 134}
135
136//________________________________________________________________________
914d486c 137void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 138{
139 // Find jets.
140
d03084cd 141 if (!tracks && !clus)
142 return;
143
48d874ff 144 TString name("kt");
145 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
146 if (algo>=1) {
147 name = "antikt";
148 jalgo = fastjet::antikt_algorithm;
149 }
150
151 AliFJWrapper fjw(name, name);
f5f3c8e9 152 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
48d874ff 153 fjw.SetR(radius);
6ea168d0 154 fjw.SetAlgorithm(jalgo);
48d874ff 155 fjw.SetMaxRap(0.9);
156 fjw.Clear();
157
158 if (tracks) {
159 const Int_t Ntracks = tracks->GetEntries();
160 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
7efbea04 161 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
48d874ff 162 if (!t)
163 continue;
914d486c 164 if (t->Pt()<fMinJetTrackPt)
165 continue;
48d874ff 166 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
167 }
168 }
6ea168d0 169
48d874ff 170 if (clus) {
171 Double_t vertex[3] = {0, 0, 0};
172 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
173 const Int_t Nclus = clus->GetEntries();
914d486c 174 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 175 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
d03084cd 176 if (!c)
177 continue;
48d874ff 178 if (!c->IsEMCAL())
179 continue;
180 TLorentzVector nPart;
181 c->GetMomentum(nPart, vertex);
63206072 182 Double_t et = nPart.Pt();
183 if (et<fMinJetClusPt)
914d486c 184 continue;
f5f3c8e9 185 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-123);
48d874ff 186 }
187 }
188
7efbea04 189 // run jet finder
48d874ff 190 fjw.Run();
f5f3c8e9 191
192 // get geometry
193 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
194 if (!geom) {
195 AliFatal("Can not create geometry");
196 return;
197 }
6ea168d0 198
48d874ff 199 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
200 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
f5f3c8e9 201 if (jets_incl[ij].perp()<fMinJetPt)
202 continue;
203 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 204 continue;
914d486c 205 AliEmcalJet *jet = new ((*fJets)[jetCount])
206 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
6ea168d0 207 Double_t vertex[3] = {0, 0, 0};
208 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 209 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
35789a2d 210 Int_t nt = 0;
211 Int_t nc = 0;
f5f3c8e9 212 Double_t neutralE = 0;
213 Double_t maxTrack = 0;
214 Double_t maxCluster = 0;
215 Int_t gall = 0;
216 Int_t gemc = 0;
217 jet->SetNumberOfTracks(constituents.size());
218 jet->SetNumberOfClusters(constituents.size());
48d874ff 219 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
220 Int_t uid = constituents[ic].user_index();
f5f3c8e9 221 if ((uid==-1) && (constituents[ic].kt2()<1e-25)) { //ghost particle
222 ++gall;
223 Double_t gphi = constituents[ic].phi()*TMath::RadToDeg();
224 Double_t geta = constituents[ic].eta();
225 if ((gphi>geom->GetArm1PhiMin()) && (gphi<geom->GetArm1PhiMax()) &&
226 (geta>geom->GetArm1EtaMin()) && (geta<geom->GetArm1EtaMax()))
227 ++gemc;
228 continue;
229 }
6ea168d0 230 if (uid>=0){
231 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
d03084cd 232 if (t) {
233 if (t->Pt()>maxTrack)
234 maxTrack=t->Pt();
f5f3c8e9 235 jet->AddTrackAt(uid, nt);
d03084cd 236 nt++;
237 }
6ea168d0 238 } else {
f5f3c8e9 239 Int_t cid = -uid - 123;
240 AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
d03084cd 241 if (c) {
242 TLorentzVector nP;
243 c->GetMomentum(nP, vertex);
244 neutralE += nP.P();
245 if (nP.Pt()>maxCluster)
246 maxCluster=nP.Pt();
f5f3c8e9 247 jet->AddClusterAt(cid,nc);
248 nc++;
d03084cd 249 }
48d874ff 250 }
251 }
35789a2d 252 jet->SetNumberOfTracks(nt);
253 jet->SetNumberOfClusters(nc);
f5f3c8e9 254 jet->SortConstituents();
6ea168d0 255 jet->SetMaxTrackPt(maxTrack);
256 jet->SetMaxClusterPt(maxCluster);
48d874ff 257 jet->SetNEF(neutralE/jet->E());
f5f3c8e9 258 jet->SetArea(fjw.GetJetArea(ij));
259 if (gall>0)
260 jet->SetAreaEmc(fjw.GetJetArea(ij)*gemc/gall);
261 else
262 jet->SetAreaEmc(-1);
263 if ((jet->Phi()>geom->GetArm1PhiMin()*TMath::DegToRad()) && (jet->Phi()<geom->GetArm1PhiMax()*TMath::DegToRad()) &&
264 (jet->Eta()>geom->GetArm1EtaMin()) && (jet->Eta()<geom->GetArm1EtaMax()))
265 jet->SetAxisInEmcal(kTRUE);
48d874ff 266 jetCount++;
267 }
268}