e6a362f6e61b041d079767109186c1657dc6a562
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx
1 // $Id$
2 //
3 // Emcal jet finder task.
4 //
5 // Authors: C.Loizides, S.Aiola
6
7 #include <vector>
8 #include "AliEmcalJetTask.h"
9
10 #include <TChain.h>
11 #include <TClonesArray.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14
15 #include "AliAnalysisManager.h"
16 #include "AliCentrality.h"
17 #include "AliEMCALGeometry.h"
18 #include "AliESDEvent.h"
19 #include "AliEmcalJet.h"
20 #include "AliEmcalParticle.h"
21 #include "AliFJWrapper.h"
22 #include "AliMCEvent.h"
23 #include "AliVCluster.h"
24 #include "AliVEvent.h"
25 #include "AliVParticle.h"
26
27 ClassImp(AliEmcalJetTask)
28
29 //________________________________________________________________________
30 AliEmcalJetTask::AliEmcalJetTask() : 
31   AliAnalysisTaskSE("AliEmcalJetTask"),
32   fTracksName("Tracks"),
33   fCaloName("CaloClusters"),
34   fJetsName("Jets"),
35   fAlgo(1),
36   fRadius(0.4),
37   fType(0),
38   fMinJetTrackPt(0.15),
39   fMinJetClusPt(0.15),
40   fPhiMin(0),
41   fPhiMax(TMath::TwoPi()),
42   fEtaMin(-1),
43   fEtaMax(1),
44   fMinJetArea(0.01),
45   fMinJetPt(1.0),
46   fIsInit(0),
47   fIsMcPart(0),
48   fIsEmcPart(0),
49   fJets(0),
50   fEvent(0),
51   fTracks(0),
52   fClus(0)
53 {
54   // Default constructor.
55 }
56
57 //________________________________________________________________________
58 AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
59   AliAnalysisTaskSE(name),
60   fTracksName("Tracks"),
61   fCaloName("CaloClusters"),
62   fJetsName("Jets"),
63   fAlgo(1),
64   fRadius(0.4),
65   fType(0),
66   fMinJetTrackPt(0.15),
67   fMinJetClusPt(0.15),
68   fPhiMin(0),
69   fPhiMax(TMath::TwoPi()),
70   fEtaMin(-1),
71   fEtaMax(1),
72   fMinJetArea(0.01),
73   fMinJetPt(1.0),
74   fIsInit(0),
75   fIsMcPart(0),
76   fIsEmcPart(0),
77   fJets(0),
78   fEvent(0),
79   fTracks(0),
80   fClus(0)
81 {
82   // Standard constructor.
83
84   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
85 }
86
87 //________________________________________________________________________
88 AliEmcalJetTask::~AliEmcalJetTask()
89 {
90   // Destructor
91 }
92
93 //________________________________________________________________________
94 void AliEmcalJetTask::UserCreateOutputObjects()
95 {
96   // Create user objects.
97
98   fJets = new TClonesArray("AliEmcalJet");
99   fJets->SetName(fJetsName);
100 }
101
102 //________________________________________________________________________
103 void AliEmcalJetTask::UserExec(Option_t *) 
104 {
105   // Main loop, called for each event.
106
107   if (!fIsInit) {
108     if (!DoInit())
109       return;
110     fIsInit = kTRUE;
111   }
112
113   FindJets();
114 }
115
116 //________________________________________________________________________
117 void AliEmcalJetTask::Terminate(Option_t *) 
118 {
119   // Called once at the end of the analysis.
120 }
121
122 //________________________________________________________________________
123 void AliEmcalJetTask::FindJets()
124 {
125   // Find jets.
126
127   if (!fTracks && !fClus)
128     return;
129
130   TString name("kt");
131   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
132   if (fAlgo>=1) {
133     name  = "antikt";
134     jalgo = fastjet::antikt_algorithm;
135   }
136
137   // setup fj wrapper
138   AliFJWrapper fjw(name, name);
139   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
140   fjw.SetR(fRadius);
141   fjw.SetAlgorithm(jalgo);  
142   fjw.SetMaxRap(1);
143   fjw.Clear();
144
145   // get primary vertex
146   Double_t vertex[3] = {0, 0, 0};
147   fEvent->GetPrimaryVertex()->GetXYZ(vertex);
148
149   if ((fIsMcPart || fType == 0 || fType == 1) && fTracks) {
150     const Int_t Ntracks = fTracks->GetEntries();
151     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
152       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
153       if (!t)
154         continue;
155       if (fIsMcPart && fType == 1 && t->Charge() == 0)
156         continue;
157       if (fIsMcPart && fType == 2 && t->Charge() != 0)
158         continue;
159       if (t->Pt() < fMinJetTrackPt) 
160         continue;
161       Double_t eta = t->Eta();
162       Double_t phi = t->Phi();
163       if ((eta<fEtaMin) && (eta>fEtaMax) &&
164           (phi<fPhiMin) && (phi<fPhiMax))
165         continue;
166
167       // offset of 100 for consistency with cluster ids
168       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  
169     }
170   }
171
172   if ((fType == 0 || fType == 2) && fClus) {
173     const Int_t Nclus = fClus->GetEntries();
174     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
175       AliVCluster *c = 0;
176       Double_t cEta=0,cPhi=0,cPt=0;
177       Double_t cPx=0,cPy=0,cPz=0;
178       if (fIsEmcPart) {
179         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
180         if (!ep)
181           continue;
182         c = ep->GetCluster();
183         if (!c)
184             continue;
185         cEta = ep->Eta();
186         cPhi = ep->Phi();
187         cPt  = ep->Pt();
188         cPx  = ep->Px();
189         cPy  = ep->Py();
190         cPz  = ep->Pz();
191       } else {
192         c = static_cast<AliVCluster*>(fClus->At(iClus));
193         if (!c)
194           continue;
195         TLorentzVector nP;
196         c->GetMomentum(nP, vertex);
197         cEta = nP.Eta();
198         cPhi = nP.Phi();
199         cPt  = nP.Pt();
200         cPx  = nP.Px();
201         cPy  = nP.Py();
202         cPz  = nP.Pz();
203       }
204       if (!c->IsEMCAL())
205         continue;
206       if (cPt < fMinJetClusPt) 
207         continue;
208       if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
209           (cPhi<fPhiMin) && (cPhi<fPhiMax))
210         continue;
211       // offset of 100 to skip ghost particles uid = -1
212       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);   
213     }
214   }
215   
216   // run jet finder
217   fjw.Run();
218
219   // get geometry
220   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
221   if (!geom) {
222     AliFatal(Form("%s: Can not create geometry", GetName()));
223     return;
224   }
225
226   // loop over fastjet jets
227   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
228   for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
229     if (jets_incl[ij].perp()<fMinJetPt) 
230       continue;
231     if (fjw.GetJetArea(ij)<fMinJetArea)
232       continue;
233
234     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
235       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
236
237     // loop over constituents
238     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
239     jet->SetNumberOfTracks(constituents.size());
240     jet->SetNumberOfClusters(constituents.size());
241
242     Int_t nt            = 0;
243     Int_t nc            = 0;
244     Double_t neutralE   = 0;
245     Double_t maxCh      = 0;
246     Double_t maxNe      = 0;
247     Int_t gall          = 0;
248     Int_t gemc          = 0;
249     Int_t cemc          = 0;
250     Int_t ncharged      = 0;
251     Int_t nneutral      = 0;
252     Double_t mcpt       = 0;
253     Double_t emcpt      = 0;
254
255     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
256       Int_t uid = constituents[ic].user_index();
257
258       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
259         ++gall;
260         Double_t gphi = constituents[ic].phi();
261         if (gphi<0) 
262           gphi += TMath::TwoPi();
263         gphi *= TMath::RadToDeg();
264         Double_t geta = constituents[ic].eta();
265         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
266             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
267           ++gemc; 
268       } else if ((uid > 0) && fTracks) { // track constituent
269         Int_t tid = uid - 100;
270         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
271         if (!t)
272           continue;
273         Double_t cEta = t->Eta();
274         Double_t cPhi = t->Phi();
275         Double_t cPt  = t->Pt();
276         Double_t cP   = t->P();
277         if (t->Charge() == 0) {
278           neutralE += cP;
279           ++nneutral;
280           if (cPt > maxNe)
281             maxNe = cPt;
282         } else {
283           ++ncharged;
284           if (cPt > maxCh)
285             maxCh = cPt;
286         }
287
288         if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
289           mcpt += cPt;
290
291         if (cPhi<0) 
292           cPhi += TMath::TwoPi();
293         cPhi *= TMath::RadToDeg();
294         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
295             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
296           emcpt += cPt;
297           ++cemc;
298         }
299
300         jet->AddTrackAt(tid, nt);
301         ++nt;
302       } else if (fClus) { // cluster constituent
303         Int_t cid = -uid - 100;
304         AliVCluster *c = 0;
305         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
306         if (fIsEmcPart) {
307           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
308           if (!ep)
309             continue;
310           c = ep->GetCluster();
311           if (!c)
312             continue;
313           cEta = ep->Eta();
314           cPhi = ep->Phi();
315           cPt  = ep->Pt();
316           cP   = ep->P();
317         } else {
318           c = static_cast<AliVCluster*>(fClus->At(cid));
319           if (!c)
320             continue;
321           TLorentzVector nP;
322           c->GetMomentum(nP, vertex);
323           cEta = nP.Eta();
324           cPhi = nP.Phi();
325           cPt  = nP.Pt();
326           cP   = nP.P();
327         }
328
329         neutralE += cP;
330         if (cPt > maxNe)
331           maxNe = cPt;
332
333         if (c->Chi2() == 100) // MC particle
334           mcpt += cPt;
335
336         if (cPhi<0) 
337           cPhi += TMath::TwoPi();
338         cPhi *= TMath::RadToDeg();
339         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
340             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
341           emcpt += cPt;
342           ++cemc;
343         }
344
345         jet->AddClusterAt(cid, nc);
346         ++nc;
347         ++nneutral;
348       } else {
349         AliError(Form("%s: No logical way to end up here.", GetName()));
350         continue;
351       }
352     } /* end of constituent loop */
353
354     jet->SetNumberOfTracks(nt);
355     jet->SetNumberOfClusters(nc);
356     jet->SortConstituents();
357     jet->SetMaxNeutralPt(maxNe);
358     jet->SetMaxChargedPt(maxCh);
359     jet->SetNEF(neutralE / jet->E());
360     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
361     jet->SetArea(area.perp());
362     jet->SetAreaEta(area.eta());
363     jet->SetAreaPhi(area.phi());
364     jet->SetNumberOfCharged(ncharged);
365     jet->SetNumberOfNeutrals(nneutral);
366     jet->SetMCPt(mcpt);
367     jet->SetNEmc(cemc);
368     jet->SetPtEmc(emcpt);
369
370     if (gall > 0)
371       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
372     else 
373       jet->SetAreaEmc(-1);
374     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
375         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
376         (jet->Eta() > geom->GetArm1EtaMin()) && 
377         (jet->Eta() < geom->GetArm1EtaMax()))
378       jet->SetAxisInEmcal(kTRUE);
379     jetCount++;
380   }
381   fJets->Sort();
382 }
383
384 //________________________________________________________________________
385 Bool_t AliEmcalJetTask::DoInit()
386 {
387   // Init. Return true if successful.
388
389   // get input collections
390   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
391
392   // get the event
393   fEvent = InputEvent();
394   if (!fEvent) {
395     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
396     return 0;
397   }
398
399   // add jets to event if not yet there
400   if (!(fEvent->FindListObject(fJetsName)))
401     fEvent->AddObject(fJets);
402   else {
403     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
404     return 0;
405   }
406
407   if (fTracksName == "Tracks")
408     am->LoadBranch("Tracks");
409   if (!fTracks && !fTracksName.IsNull()) {
410     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
411     if (!fTracks) {
412       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
413       return 0;
414     }
415   }
416   if (fTracks) {
417     TClass cls(fTracks->GetClass()->GetName());
418     if (cls.InheritsFrom("AliMCParticle"))
419       fIsMcPart = 1;
420   }
421   
422   if (fCaloName == "CaloClusters")
423     am->LoadBranch("CaloClusters");
424   if (!fClus && !fCaloName.IsNull()) {
425     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
426     if (!fClus) {
427       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
428       return 0;
429     }
430   }
431   if (fClus) {
432     TClass cls(fClus->GetClass()->GetName());
433     if (cls.InheritsFrom("AliEmcalParticle"))
434       fIsEmcPart = 1;
435   }
436
437   return 1;
438 }