Compute jet area in emcal from ghost particles. Introduce handy getters for IsInsideE...
[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"
63206072 17#include "AliFJWrapper.h"
35789a2d 18#include "AliVCluster.h"
19#include "AliVTrack.h"
914d486c 20#include "AliEmcalJet.h"
48d874ff 21
914d486c 22ClassImp(AliEmcalJetTask)
48d874ff 23
48d874ff 24//________________________________________________________________________
914d486c 25AliEmcalJetTask::AliEmcalJetTask() :
26 AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0 27 fTracksName("Tracks"),
28 fCaloName("CaloClusters"),
29 fJetsName("Jets"),
30 fAlgo(1),
31 fRadius(0.4),
32 fType(0),
914d486c 33 fMinJetTrackPt(0.15),
34 fMinJetClusPt(0.15),
6ea168d0 35 fJets(0)
36{
37 // Default constructor.
6ea168d0 38}
39
40//________________________________________________________________________
914d486c 41AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072 42 AliAnalysisTaskSE(name),
7efbea04 43 fTracksName("Tracks"),
48d874ff 44 fCaloName("CaloClusters"),
7efbea04 45 fJetsName("Jets"),
48d874ff 46 fAlgo(1),
47 fRadius(0.4),
48 fType(0),
6ea168d0 49 fMinJetTrackPt(0.15),
914d486c 50 fMinJetClusPt(0.15),
7efbea04 51 fJets(0)
48d874ff 52{
53 // Standard constructor.
7efbea04 54
7efbea04 55 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff 56}
57
58//________________________________________________________________________
914d486c 59AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 60{
61 // Destructor
62}
63
64//________________________________________________________________________
914d486c 65void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff 66{
67 // Create user objects.
68
914d486c 69 fJets = new TClonesArray("AliEmcalJet");
48d874ff 70 fJets->SetName(fJetsName);
48d874ff 71}
72
73//________________________________________________________________________
914d486c 74void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 75{
76 // Main loop, called for each event.
77
63206072 78 // add jets to event if not yet there
48d874ff 79 if (!(InputEvent()->FindListObject(fJetsName)))
80 InputEvent()->AddObject(fJets);
81
63206072 82 // delete jet output
83 fJets->Delete();
84
85
86 // get input collections
48d874ff 87 TClonesArray *tracks = 0;
88 TClonesArray *clus = 0;
89 TList *l = InputEvent()->GetList();
63206072 90 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
6ea168d0 91
48d874ff 92 if ((fType==0)||(fType==1)) {
7efbea04 93 if (fTracksName == "Tracks")
48d874ff 94 am->LoadBranch("Tracks");
7efbea04 95 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
48d874ff 96 if (!tracks) {
653159b9 97 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
48d874ff 98 return;
99 }
100 }
101 if ((fType==0)||(fType==2)) {
102 if (fCaloName == "CaloClusters")
103 am->LoadBranch("CaloClusters");
104 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
105 if (!clus) {
653159b9 106 AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
48d874ff 107 return;
108 }
48d874ff 109 }
63206072 110
111 // get centrality
112 Float_t cent = -1;
653159b9 113 AliCentrality *centrality = InputEvent()->GetCentrality();
63206072 114 if (centrality)
115 cent = centrality->GetCentralityPercentile("V0M");
116 else
117 cent=99; // probably pp data
653159b9 118
63206072 119 if (cent<0) {
120 AliError(Form("Centrality negative: %f", cent));
121 return;
122 }
653159b9 123
6ea168d0 124 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 125}
126
127//________________________________________________________________________
914d486c 128void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 129{
130 // Called once at the end of the analysis.
48d874ff 131}
132
133//________________________________________________________________________
914d486c 134void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 135{
136 // Find jets.
137
d03084cd 138 if (!tracks && !clus)
139 return;
140
48d874ff 141 TString name("kt");
142 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
143 if (algo>=1) {
144 name = "antikt";
145 jalgo = fastjet::antikt_algorithm;
146 }
147
148 AliFJWrapper fjw(name, name);
149 fjw.SetR(radius);
6ea168d0 150 fjw.SetAlgorithm(jalgo);
48d874ff 151 fjw.SetMaxRap(0.9);
152 fjw.Clear();
153
154 if (tracks) {
155 const Int_t Ntracks = tracks->GetEntries();
156 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
7efbea04 157 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
48d874ff 158 if (!t)
159 continue;
914d486c 160 if (t->Pt()<fMinJetTrackPt)
161 continue;
48d874ff 162 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
163 }
164 }
6ea168d0 165
48d874ff 166 if (clus) {
167 Double_t vertex[3] = {0, 0, 0};
168 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
169 const Int_t Nclus = clus->GetEntries();
914d486c 170 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 171 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
d03084cd 172 if (!c)
173 continue;
48d874ff 174 if (!c->IsEMCAL())
175 continue;
176 TLorentzVector nPart;
177 c->GetMomentum(nPart, vertex);
63206072 178 Double_t et = nPart.Pt();
179 if (et<fMinJetClusPt)
914d486c 180 continue;
63206072 181 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-1);
48d874ff 182 }
183 }
184
7efbea04 185 // run jet finder
48d874ff 186 fjw.Run();
6ea168d0 187
48d874ff 188 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
189 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
190 if (jets_incl[ij].perp()<1)
191 continue;
914d486c 192 AliEmcalJet *jet = new ((*fJets)[jetCount])
193 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
900f5135 194 jet->SetArea(fjw.GetJetArea(ij));
6ea168d0 195 Double_t vertex[3] = {0, 0, 0};
196 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 197 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
6ea168d0 198 Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
35789a2d 199 jet->SetNumberOfTracks(constituents.size());
200 jet->SetNumberOfClusters(constituents.size());
201 Int_t nt = 0;
202 Int_t nc = 0;
48d874ff 203 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
204 Int_t uid = constituents[ic].user_index();
6ea168d0 205 if (uid>=0){
35789a2d 206 jet->AddTrackAt(uid, nt);
6ea168d0 207 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
d03084cd 208 if (t) {
209 if (t->Pt()>maxTrack)
210 maxTrack=t->Pt();
211 nt++;
212 }
6ea168d0 213 } else {
35789a2d 214 jet->AddClusterAt(-(uid+1),nc);
d03084cd 215 AliVCluster *c = static_cast<AliVCluster*>(clus->At(-(uid+1)));
216 if (c) {
217 TLorentzVector nP;
218 c->GetMomentum(nP, vertex);
219 neutralE += nP.P();
220 if (nP.Pt()>maxCluster)
221 maxCluster=nP.Pt();
222 }
35789a2d 223 nc++;
48d874ff 224 }
225 }
35789a2d 226 jet->SetNumberOfTracks(nt);
227 jet->SetNumberOfClusters(nc);
6ea168d0 228 jet->SetMaxTrackPt(maxTrack);
229 jet->SetMaxClusterPt(maxCluster);
48d874ff 230 jet->SetNEF(neutralE/jet->E());
48d874ff 231 jetCount++;
232 }
233}