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