c5f86d314deddfc0a407618aa5b61970c9f7737f
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliEmcalJetTask.cxx
1 // $Id$
2 //
3 // Emcal jet finder task.
4 //
5 // Author: C.Loizides
6
7 #include "AliEmcalJetTask.h"
8 #include <TChain.h>
9 #include <TClonesArray.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TParticle.h>
15 #include "AliAnalysisManager.h"
16 #include "AliCentrality.h"
17 #include "AliEMCALGeometry.h"
18 #include "AliEmcalJet.h"
19 #include "AliFJWrapper.h"
20 #include "AliVCluster.h"
21 #include "AliVTrack.h"
22
23 ClassImp(AliEmcalJetTask)
24
25 //________________________________________________________________________
26 AliEmcalJetTask::AliEmcalJetTask() : 
27   AliAnalysisTaskSE("AliEmcalJetTask"),
28   fTracksName("Tracks"),
29   fCaloName("CaloClusters"),
30   fJetsName("Jets"),
31   fAlgo(1),
32   fRadius(0.4),
33   fType(0),
34   fMinJetTrackPt(0.15),
35   fMinJetClusPt(0.15),
36   fMinJetArea(0.01),
37   fMinJetPt(1.0),
38   fJets(0)
39 {
40   // Default constructor.
41 }
42
43 //________________________________________________________________________
44 AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
45   AliAnalysisTaskSE(name),
46   fTracksName("Tracks"),
47   fCaloName("CaloClusters"),
48   fJetsName("Jets"),
49   fAlgo(1),
50   fRadius(0.4),
51   fType(0),
52   fMinJetTrackPt(0.15),
53   fMinJetClusPt(0.15),
54   fMinJetArea(0.01),
55   fMinJetPt(1.0),
56   fJets(0)
57 {
58   // Standard constructor.
59
60   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
61 }
62
63 //________________________________________________________________________
64 AliEmcalJetTask::~AliEmcalJetTask()
65 {
66   // Destructor
67 }
68
69 //________________________________________________________________________
70 void AliEmcalJetTask::UserCreateOutputObjects()
71 {
72   // Create user objects.
73
74   fJets = new TClonesArray("AliEmcalJet");
75   fJets->SetName(fJetsName);
76 }
77
78 //________________________________________________________________________
79 void AliEmcalJetTask::UserExec(Option_t *) 
80 {
81   // Main loop, called for each event.
82
83   // add jets to event if not yet there
84   if (!(InputEvent()->FindListObject(fJetsName)))
85     InputEvent()->AddObject(fJets);
86
87   // delete jet output
88   fJets->Delete();
89
90   // get input collections
91   TClonesArray *tracks = 0;
92   TClonesArray *clus   = 0;
93   TList *l = InputEvent()->GetList();
94   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
95
96   if ((fType==0)||(fType==1)) {
97     if (fTracksName == "Tracks")
98       am->LoadBranch("Tracks");
99     tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
100     if (!tracks) {
101       AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
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) {
110       AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
111       return;
112     }
113   }
114
115   // get centrality
116   Float_t cent = -1; 
117   AliCentrality *centrality = InputEvent()->GetCentrality();
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   }
126
127   FindJets(tracks, clus, fAlgo, fRadius, cent);
128 }
129
130 //________________________________________________________________________
131 void AliEmcalJetTask::Terminate(Option_t *) 
132 {
133   // Called once at the end of the analysis.
134 }
135
136 //________________________________________________________________________
137 void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
138 {
139   // Find jets.
140
141   if (!tracks && !clus)
142     return;
143
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);
152   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
153   fjw.SetR(radius);
154   fjw.SetAlgorithm(jalgo);  
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) {
161       AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
162       if (!t)
163         continue;
164       if (t->Pt()<fMinJetTrackPt) 
165         continue;
166       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
167     }
168   }
169
170   if (clus) {
171     Double_t vertex[3] = {0, 0, 0};
172     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
173     const Int_t Nclus = clus->GetEntries();
174     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
175       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
176       if (!c)
177         continue;
178       if (!c->IsEMCAL())
179         continue;
180       TLorentzVector nPart;
181       c->GetMomentum(nPart, vertex);
182       Double_t et = nPart.Pt();
183       if (et<fMinJetClusPt) 
184         continue;
185       fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-123);
186     }
187   }
188
189   // run jet finder
190   fjw.Run();
191
192   // get geometry
193   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
194   if (!geom) {
195     AliFatal("Can not create geometry");
196     return;
197   }
198   
199   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
200   for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
201     if (jets_incl[ij].perp()<fMinJetPt) 
202       continue;
203     if (fjw.GetJetArea(ij)<fMinJetArea)
204       continue;
205     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
206       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
207     Double_t vertex[3] = {0, 0, 0};
208     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
209     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
210     Int_t nt = 0;
211     Int_t nc = 0;
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());
219     for(UInt_t ic=0; ic<constituents.size(); ++ic) {
220       Int_t uid = constituents[ic].user_index();
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       }
230       if (uid>=0){
231         AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
232         if (t) {
233           if (t->Pt()>maxTrack)
234             maxTrack=t->Pt();
235           jet->AddTrackAt(uid, nt);
236           nt++;
237         }
238       } else {
239         Int_t cid = -uid - 123;
240         AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
241         if (c) {
242           TLorentzVector nP;
243           c->GetMomentum(nP, vertex);
244           neutralE += nP.P();
245           if (nP.Pt()>maxCluster)
246             maxCluster=nP.Pt();
247           jet->AddClusterAt(cid,nc);
248           nc++;
249         }
250       }
251     }
252     jet->SetNumberOfTracks(nt);
253     jet->SetNumberOfClusters(nc);
254     jet->SortConstituents();
255     jet->SetMaxTrackPt(maxTrack);
256     jet->SetMaxClusterPt(maxCluster);
257     jet->SetNEF(neutralE/jet->E());
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);
266     jetCount++;
267   }
268 }