have jet task and wrapper
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliEsdJetTask.cxx
CommitLineData
48d874ff 1// $Id$
2
3#include "AliEsdJetTask.h"
4#include <TClonesArray.h>
5#include <TParticle.h>
6#include "AliAODJet.h"
7#include "AliAnalysisManager.h"
8#include "AliESDtrack.h"
9#include "AliESDtrackCuts.h"
10#include "AliFJWrapper.h"
11#include "AliESDCaloCluster.h"
12
13ClassImp(AliEsdJetTask)
14
15#ifdef __HAVE_FJINTERFACE__
16//________________________________________________________________________
17AliEsdJetTask::AliEsdJetTask(const char *name) :
18 AliAnalysisTaskSE("AliEsdJetTask"),
19 fPrimTrCuts(0),
20 fPrimTracksName("Tracks"),
21 fJetsName("Jets"),
22 fCaloName("CaloClusters"),
23 fAlgo(1),
24 fRadius(0.4),
25 fType(0),
26 fHadCorr(0),
27 fJets(0),
28 fNeutrals(0)
29{
30 // Standard constructor.
31 if (!name)
32 return;
33 SetName(name);
34 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,CaloClusters,Tracks";
35}
36
37//________________________________________________________________________
38AliEsdJetTask::~AliEsdJetTask()
39{
40 // Destructor
41}
42
43//________________________________________________________________________
44void AliEsdJetTask::UserCreateOutputObjects()
45{
46 // Create user objects.
47
48 fJets = new TClonesArray("AliAODJet");
49 fJets->SetName(fJetsName);
50
51 if ((fType==0)||(fType==2)) {
52 fNeutrals = new TClonesArray("TLorentzVector");
53 fNeutrals->SetName(Form("%s_neutrals",fJetsName.Data()));
54 }
55}
56
57//________________________________________________________________________
58void AliEsdJetTask::UserExec(Option_t *)
59{
60 // Main loop, called for each event.
61
62 // add jets to event if not yet there
63 if (!(InputEvent()->FindListObject(fJetsName)))
64 InputEvent()->AddObject(fJets);
65
66 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
67 TClonesArray *tracks = 0;
68 TClonesArray *clus = 0;
69 TList *l = InputEvent()->GetList();
70 if ((fType==0)||(fType==1)) {
71 if (fPrimTracksName == "Tracks")
72 am->LoadBranch("Tracks");
73 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fPrimTracksName));
74 if (!tracks) {
75 AliError(Form("Pointer to tracks %s == 0", fPrimTracksName.Data() ));
76 return;
77 }
78 }
79 if ((fType==0)||(fType==2)) {
80 if (fCaloName == "CaloClusters")
81 am->LoadBranch("CaloClusters");
82 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
83 if (!clus) {
84 AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
85 return;
86 }
87 if (!(InputEvent()->FindListObject(Form("%s_neutrals",fJetsName.Data()))))
88 InputEvent()->AddObject(fNeutrals);
89 }
90
91 FindJets(tracks, clus, fAlgo, fRadius);
92}
93
94//________________________________________________________________________
95void AliEsdJetTask::Terminate(Option_t *)
96{
97 // Called once at the end of the analysis.
98
99}
100
101//________________________________________________________________________
102void AliEsdJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius)
103{
104 // Find jets.
105
106 TString name("kt");
107 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
108 if (algo>=1) {
109 name = "antikt";
110 jalgo = fastjet::antikt_algorithm;
111 }
112
113 AliFJWrapper fjw(name, name);
114 fjw.SetR(radius);
115 fjw.SetAlgorithm(jalgo);
116 fjw.SetMaxRap(0.9);
117 fjw.Clear();
118
119 if (tracks) {
120 const Int_t Ntracks = tracks->GetEntries();
121 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
122 AliESDtrack *t = static_cast<AliESDtrack*>(tracks->At(iTracks));
123 if (!t)
124 continue;
125 if (fPrimTrCuts && !fPrimTrCuts->IsSelected(t))
126 continue;
127 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
128 }
129 }
130 if (clus) {
131 Double_t vertex[3] = {0, 0, 0};
132 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
133 const Int_t Nclus = clus->GetEntries();
134 for (Int_t iClus = 0, iN = 0; iClus < Nclus; ++iClus) {
135 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(iClus));
136 if (!c->IsEMCAL())
137 continue;
138 TLorentzVector nPart;
139 c->GetMomentum(nPart, vertex);
140 Double_t energy = nPart.P();
141 if (fHadCorr>0) {
142 Int_t imin = static_cast<Int_t>(c->GetEmcCpvDistance());
143 if (imin>=0) {
144 Double_t dPhiMin = c->GetTrackDx();
145 Double_t dEtaMin = c->GetTrackDz();
146 if (dPhiMin<0.05 && dEtaMin<0.025) {
147 AliESDtrack *t = static_cast<AliESDtrack*>(tracks->At(imin));
148 if (t) {
149 if (!fPrimTrCuts || fPrimTrCuts->IsSelected(t)) {
150 energy -= fHadCorr*t->P();
151 if (energy<0)
152 continue;
153 }
154 }
155 }
156 }
157 }
158 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), energy, -iN-1);
159 new ((*fNeutrals)[iN++]) TLorentzVector(nPart);
160 }
161 }
162
163 fjw.Run();
164 if (0) {
165 Double_t median = 0;
166 Double_t sigma = 0;
167 fjw.GetMedianAndSigma(median, sigma);
168 }
169
170 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
171 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
172 if (jets_incl[ij].perp()<1)
173 continue;
174 AliAODJet *jet = new ((*fJets)[jetCount])
175 AliAODJet(jets_incl[ij].px(), jets_incl[ij].py(), jets_incl[ij].pz(), jets_incl[ij].E());
176 jet->SetEffArea(fjw.GetJetArea(ij),0);
177 jet->GetRefTracks()->Clear();
178 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
179 Double_t neutralE = 0;
180 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
181 Int_t uid = constituents[ic].user_index();
182 if (uid>=0)
183 jet->AddTrack(tracks->At(uid));
184 else {
185 TLorentzVector *nP = static_cast<TLorentzVector*>(fNeutrals->At(-(uid+1)));
186 neutralE += nP->E();
187 jet->AddTrack(nP);
188 }
189 }
190 jet->SetNEF(neutralE/jet->E());
191 //jet->Print("");
192 jetCount++;
193 }
194}
195#else
196//________________________________________________________________________
197AliEsdJetTask::AliEsdJetTask(const char *name) :
198 AliAnalysisTaskSE("AliEsdJetTask"),
199 fPrimTrCuts(0),
200 fPrimTracksName("Tracks"),
201 fJetsName("Jets"),
202 fCaloName("CaloClusters"),
203 fAlgo(1),
204 fRadius(0.4),
205 fType(0),
206 fHadCorr(0),
207 fJets(0),
208 fNeutrals(0)
209{
210 // Standard constructor.
211 AliFatal("Compiled without FASTJET package. Task cannot be used!");
212}
213
214//________________________________________________________________________
215AliEsdJetTask::~AliEsdJetTask()
216{
217 // Destructor
218}
219
220//________________________________________________________________________
221void AliEsdJetTask::UserCreateOutputObjects()
222{
223 // Create user objects.
224}
225
226//________________________________________________________________________
227void AliEsdJetTask::UserExec(Option_t *)
228{
229}
230
231//________________________________________________________________________
232void AliEsdJetTask::Terminate(Option_t *)
233{
234 // Called once at the end of the analysis.
235}
236
237//________________________________________________________________________
238void AliEsdJetTask::FindJets(TObjArray *, TObjArray *, Int_t, Double_t)
239{
240 // Find jets.
241}
242#endif