Compatibility with the Root trunk
[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 <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   // delete jet output
114   fJets->Delete();
115
116   FindJets();
117 }
118
119 //________________________________________________________________________
120 void AliEmcalJetTask::Terminate(Option_t *) 
121 {
122   // Called once at the end of the analysis.
123 }
124
125 //________________________________________________________________________
126 void AliEmcalJetTask::FindJets()
127 {
128   // Find jets.
129
130   if (!fTracks && !fClus)
131     return;
132
133   TString name("kt");
134   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
135   if (fAlgo>=1) {
136     name  = "antikt";
137     jalgo = fastjet::antikt_algorithm;
138   }
139
140   // setup fj wrapper
141   AliFJWrapper fjw(name, name);
142   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
143   fjw.SetR(fRadius);
144   fjw.SetAlgorithm(jalgo);  
145   fjw.SetMaxRap(1);
146   fjw.Clear();
147
148   // get primary vertex
149   Double_t vertex[3] = {0, 0, 0};
150   fEvent->GetPrimaryVertex()->GetXYZ(vertex);
151
152   if (fTracks) {
153     const Int_t Ntracks = fTracks->GetEntries();
154     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
155       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
156       if (!t)
157         continue;
158       if (t->Pt() < fMinJetTrackPt) 
159         continue;
160       Double_t eta = t->Eta();
161       Double_t phi = t->Phi();
162       if ((eta<fEtaMin) && (eta>fEtaMax) &&
163           (phi<fPhiMin) && (phi<fPhiMax))
164         continue;
165       // offset of 100 for consistency with cluster ids
166       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  
167     }
168   }
169
170   if (fClus) {
171     if (fIsEmcPart) {
172       const Int_t Nclus = fClus->GetEntries();
173       for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
174         AliVCluster *c = 0;
175         Double_t cEta=0,cPhi=0,cPt=0;
176         Double_t cPx=0,cPy=0,cPz=0;
177         if (fIsEmcPart) {
178           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
179           if (!ep)
180             continue;
181           c = ep->GetCluster();
182           if (!c)
183             continue;
184           cEta = ep->Eta();
185           cPhi = ep->Phi();
186           cPt  = ep->Pt();
187           cPx  = ep->Px();
188           cPy  = ep->Py();
189           cPz  = ep->Pz();
190         } else {
191           c = static_cast<AliVCluster*>(fClus->At(iClus));
192           if (!c)
193             continue;
194           TLorentzVector nP;
195           c->GetMomentum(nP, vertex);
196           cEta = nP.Eta();
197           cPhi = nP.Phi();
198           cPt  = nP.Pt();
199           cPx  = nP.Px();
200           cPy  = nP.Py();
201           cPz  = nP.Pz();
202         }
203         if (!c->IsEMCAL())
204           continue;
205         if (cPt < fMinJetClusPt) 
206           continue;
207         if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
208             (cPhi<fPhiMin) && (cPhi<fPhiMax))
209           continue;
210         // offset of 100 to skip ghost particles uid = -1
211         fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);  
212       }
213     } else { /*CaloClusters given as input*/
214       const Int_t Nclus = fClus->GetEntries();
215       for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
216         AliVCluster *c = static_cast<AliVCluster*>(fClus->At(iClus));
217         if (!c)
218           continue;
219         if (!c->IsEMCAL())
220           continue;
221         TLorentzVector nPart;
222         c->GetMomentum(nPart, vertex);
223         if (nPart.Pt() < fMinJetClusPt) 
224           continue;
225         // offset of 100 to skip ghost particles uid = -1
226         fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);  
227       }
228     }
229   }
230
231   // run jet finder
232   fjw.Run();
233
234   // get geometry
235   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
236   if (!geom) {
237     AliFatal(Form("%s: Can not create geometry", GetName()));
238     return;
239   }
240
241   // loop over fastjet jets
242   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
243   for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
244     if (jets_incl[ij].perp()<fMinJetPt) 
245       continue;
246     if (fjw.GetJetArea(ij)<fMinJetArea)
247       continue;
248
249     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
250       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
251
252     // loop over constituents
253     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
254     jet->SetNumberOfTracks(constituents.size());
255     jet->SetNumberOfClusters(constituents.size());
256
257     Int_t nt            = 0;
258     Int_t nc            = 0;
259     Double_t neutralE   = 0;
260     Double_t maxCh      = 0;
261     Double_t maxNe      = 0;
262     Int_t gall          = 0;
263     Int_t gemc          = 0;
264     Int_t cemc          = 0;
265     Int_t ncharged      = 0;
266     Int_t nneutral      = 0;
267     Double_t mcpt       = 0;
268     Double_t emcpt      = 0;
269
270     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
271       Int_t uid = constituents[ic].user_index();
272
273       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
274         ++gall;
275         Double_t gphi = constituents[ic].phi();
276         if (gphi<0) 
277           gphi += TMath::TwoPi();
278         gphi *= TMath::RadToDeg();
279         Double_t geta = constituents[ic].eta();
280         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
281             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
282           ++gemc; 
283       } else if ((uid > 0) && fTracks) { // track constituent
284         Int_t tid = uid - 100;
285         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
286         if (!t)
287           continue;
288         Double_t cEta = t->Eta();
289         Double_t cPhi = t->Phi();
290         Double_t cPt  = t->Pt();
291         Double_t cP   = t->P();
292         if (t->Charge() == 0) {
293           neutralE += cP;
294           ++nneutral;
295           if (cPt > maxNe)
296             maxNe = cPt;
297         } else {
298           ++ncharged;
299           if (cPt > maxCh)
300             maxCh = cPt;
301         }
302
303         if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
304           mcpt += cPt;
305
306         if (cPhi<0) 
307           cPhi += TMath::TwoPi();
308         cPhi *= TMath::RadToDeg();
309         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
310             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
311           emcpt += cPt;
312           ++cemc;
313         }
314
315         jet->AddTrackAt(tid, nt);
316         ++nt;
317       } else if (fClus) { // cluster constituent
318         Int_t cid = -uid - 100;
319         AliVCluster *c = 0;
320         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
321         if (fIsEmcPart) {
322           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
323           if (!ep)
324             continue;
325           c = ep->GetCluster();
326           if (!c)
327             continue;
328           cEta = ep->Eta();
329           cPhi = ep->Phi();
330           cPt  = ep->Pt();
331           cP   = ep->P();
332         } else {
333           c = static_cast<AliVCluster*>(fClus->At(cid));
334           if (!c)
335             continue;
336           TLorentzVector nP;
337           c->GetMomentum(nP, vertex);
338           cEta = nP.Eta();
339           cPhi = nP.Phi();
340           cPt  = nP.Pt();
341           cP   = nP.P();
342         }
343
344         neutralE += cP;
345         if (cPt > maxNe)
346           maxNe = cPt;
347
348         if (c->Chi2() == 100) // MC particle
349           mcpt += cPt;
350
351         if (cPhi<0) 
352           cPhi += TMath::TwoPi();
353         cPhi *= TMath::RadToDeg();
354         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
355             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
356           emcpt += cPt;
357           ++cemc;
358         }
359
360         jet->AddClusterAt(cid, nc);
361         ++nc;
362         ++nneutral;
363       } else {
364         AliError(Form("%s: No logical way to end up here.", GetName()));
365         continue;
366       }
367     } /* end of constituent loop */
368
369     jet->SetNumberOfTracks(nt);
370     jet->SetNumberOfClusters(nc);
371     jet->SortConstituents();
372     jet->SetMaxNeutralPt(maxNe);
373     jet->SetMaxChargedPt(maxCh);
374     jet->SetNEF(neutralE / jet->E());
375     jet->SetArea(fjw.GetJetArea(ij));
376     jet->SetNumberOfCharged(ncharged);
377     jet->SetNumberOfNeutrals(nneutral);
378     jet->SetMCPt(mcpt);
379     jet->SetNEmc(cemc);
380     jet->SetPtEmc(emcpt);
381
382     if (gall > 0)
383       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
384     else 
385       jet->SetAreaEmc(-1);
386     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
387         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
388         (jet->Eta() > geom->GetArm1EtaMin()) && 
389         (jet->Eta() < geom->GetArm1EtaMax()))
390       jet->SetAxisInEmcal(kTRUE);
391     jetCount++;
392   }
393 }
394
395 //________________________________________________________________________
396 Bool_t AliEmcalJetTask::DoInit()
397 {
398   // Init. Return true if successful.
399
400   // get input collections
401   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
402
403   // get the event
404   fEvent = InputEvent();
405   if (!fEvent) {
406     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
407     return 0;
408   }
409
410   // add jets to event if not yet there
411   if (!(fEvent->FindListObject(fJetsName)))
412     fEvent->AddObject(fJets);
413   else {
414     AliError(Form("%s: Object with name %s already in event! Returning", fJetsName.Data(), GetName()));
415     return 0;
416   }
417
418   if ((fType == 0) || (fType == 1)) {
419     if (fTracksName == "Tracks")
420       am->LoadBranch("Tracks");
421     if (!fTracks && !fTracksName.IsNull()) {
422       fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
423       if (!fTracks) {
424         AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
425         return 0;
426       }
427     }
428     TString objname(fTracks->GetClass()->GetName());
429     TClass cls(objname);
430     if (cls.InheritsFrom("AliMCParticle"))
431       fIsMcPart = 1;
432   }
433
434   if ((fType == 0) || (fType == 2)) {
435     if (fCaloName == "CaloClusters")
436       am->LoadBranch("CaloClusters");
437     if (!fClus && !fCaloName.IsNull()) {
438       fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
439       if (!fClus) {
440         AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
441         return 0;
442       }
443     }
444     TString objname(fClus->GetClass()->GetName());
445     TClass cls(objname);
446     if (cls.InheritsFrom("AliEmcalParticle"))
447       fIsEmcPart = 1;
448   }
449
450   return 1;
451 }