]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/EMCALJetTasks/AliEmcalJetTask.cxx
coverty
[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 "AliESDEvent.h"
21 #include "AliMCEvent.h"
22
23 #include "AliEmcalJetTask.h"
24
25 ClassImp(AliEmcalJetTask)
26
27 //________________________________________________________________________
28 AliEmcalJetTask::AliEmcalJetTask() : 
29   AliAnalysisTaskSE("AliEmcalJetTask"),
30   fTracksName("Tracks"),
31   fCaloName("CaloClusters"),
32   fJetsName("Jets"),
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   fAlgo(1),
53   fRadius(0.4),
54   fType(0),
55   fMinJetTrackPt(0.15),
56   fMinJetClusPt(0.15),
57   fMinJetArea(0.01),
58   fMinJetPt(1.0),
59   fJets(0),
60   fEvent(0)
61 {
62   // Standard constructor.
63
64   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
65 }
66
67 //________________________________________________________________________
68 AliEmcalJetTask::~AliEmcalJetTask()
69 {
70   // Destructor
71 }
72
73 //________________________________________________________________________
74 void AliEmcalJetTask::UserCreateOutputObjects()
75 {
76   // Create user objects.
77
78   fJets = new TClonesArray("AliEmcalJet");
79   fJets->SetName(fJetsName);
80 }
81
82 //_____________________________________________________
83 TString AliEmcalJetTask::GetBeamType()
84 {
85   // Get beam type : pp-AA-pA
86   // ESDs have it directly, AODs get it from hardcoded run number ranges
87
88   TString beamType;
89
90   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
91   if (esd) {
92     const AliESDRun *run = esd->GetESDRun();
93     beamType = run->GetBeamType();
94   }
95   else
96   {
97     Int_t runNumber = fEvent->GetRunNumber();
98     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
99         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
100     {
101       beamType = "A-A";
102     }
103     else 
104     {
105       beamType = "p-p";
106     }
107   }
108
109   return beamType;    
110 }
111
112
113 //________________________________________________________________________
114 void AliEmcalJetTask::UserExec(Option_t *) 
115 {
116   // Main loop, called for each event.
117
118   // get the event
119   fEvent = InputEvent();
120
121   if (!fEvent) {
122     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
123     return;
124   }
125
126   // add jets to event if not yet there
127   if (!(fEvent->FindListObject(fJetsName)))
128     fEvent->AddObject(fJets);
129
130   // delete jet output
131   fJets->Delete();
132
133   // get input collections
134   TClonesArray *tracks = 0;
135   TClonesArray *clus   = 0;
136   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
137
138   if ((fType == 0) || (fType == 1)) {
139     if (fTracksName == "Tracks")
140       am->LoadBranch("Tracks");
141     tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
142     if (!tracks) {
143       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
144       return;
145     }
146   }
147   if ((fType == 0) || (fType == 2)) {
148     if (fCaloName == "CaloClusters")
149       am->LoadBranch("CaloClusters");
150     clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
151     if (!clus) {
152       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
153       return;
154     }
155   }
156   
157   // get centrality 
158   Double_t cent = 99; 
159   
160   if (GetBeamType() == "A-A") {
161     AliCentrality *centrality = InputEvent()->GetCentrality();
162     
163     if (centrality)
164       cent = centrality->GetCentralityPercentile("V0M");
165     else
166       cent = 99; // probably pp data
167     
168     if (cent < 0) {
169       AliWarning(Form("Centrality negative: %f, assuming 99", cent));
170       cent = 99;
171     }
172   }
173
174   FindJets(tracks, clus, fAlgo, fRadius, cent);
175 }
176
177 //________________________________________________________________________
178 void AliEmcalJetTask::Terminate(Option_t *) 
179 {
180   // Called once at the end of the analysis.
181 }
182
183 //________________________________________________________________________
184 void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
185 {
186   // Find jets.
187
188   if (!tracks && !clus)
189     return;
190
191   TString name("kt");
192   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
193   if (algo>=1) {
194     name  = "antikt";
195     jalgo = fastjet::antikt_algorithm;
196   }
197
198   AliFJWrapper fjw(name, name);
199   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
200   fjw.SetR(radius);
201   fjw.SetAlgorithm(jalgo);  
202   fjw.SetMaxRap(0.9);
203   fjw.Clear();
204
205   if (tracks) {
206     const Int_t Ntracks = tracks->GetEntries();
207     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
208       AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
209       if (!t)
210         continue;
211
212       if (t->Pt() < fMinJetTrackPt) 
213         continue;
214       // offset of 100 for consistency with cluster ids
215       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  
216     }
217   }
218
219   if (clus) {
220     Double_t vertex[3] = {0, 0, 0};
221     fEvent->GetPrimaryVertex()->GetXYZ(vertex);
222     const Int_t Nclus = clus->GetEntries();
223     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
224       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
225       if (!c)
226         continue;
227       if (!c->IsEMCAL())
228         continue;
229       TLorentzVector nPart;
230       c->GetMomentum(nPart, vertex);
231       if (nPart.Pt() < fMinJetClusPt) 
232         continue;
233       // offset of 100 to skip ghost particles uid = -1
234       fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);  
235     }
236   }
237
238   // run jet finder
239   fjw.Run();
240
241   // get geometry
242   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
243   if (!geom) {
244     AliFatal(Form("%s: Can not create geometry", GetName()));
245     return;
246   }
247   
248   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
249   for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
250     if (jets_incl[ij].perp()<fMinJetPt) 
251       continue;
252     if (fjw.GetJetArea(ij)<fMinJetArea)
253       continue;
254     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
255       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
256     Double_t vertex[3] = {0, 0, 0};
257     fEvent->GetPrimaryVertex()->GetXYZ(vertex);
258     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
259     Int_t nt            = 0;
260     Int_t nc            = 0;
261     Double_t neutralE   = 0;
262     Double_t maxCh      = 0;
263     Double_t maxNe      = 0;
264     Int_t gall          = 0;
265     Int_t gemc          = 0;
266     Int_t ncharged      = 0;
267     Int_t nneutral      = 0;
268     Double_t mcpt       = 0;
269
270     jet->SetNumberOfTracks(constituents.size());
271     jet->SetNumberOfClusters(constituents.size());
272     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
273       Int_t uid = constituents[ic].user_index();
274
275       if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
276         ++gall;
277         Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
278         Double_t geta = constituents[ic].eta();
279         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
280             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
281           ++gemc;
282         continue;
283       } else if ((uid > 0) && tracks) { // track constituent
284         Int_t tid = uid - 100;
285         AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
286         if (t) {
287           if (t->Charge() == 0) {
288             neutralE += t->P();
289             ++nneutral;
290             if (t->Pt() > maxNe)
291               maxNe = t->Pt();
292           } else {
293             ++ncharged;
294             if (t->Pt() > maxCh)
295               maxCh = t->Pt();
296           }
297
298           if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
299             mcpt += t->Pt();
300
301           jet->AddTrackAt(tid, nt);
302           ++nt;
303         }
304       } else if (clus) { // cluster constituent
305         Int_t cid = -uid - 100;
306         AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
307         if (c) {
308           TLorentzVector nP;
309           c->GetMomentum(nP, vertex);
310           neutralE += nP.P();
311           if (nP.Pt() > maxNe)
312             maxNe = nP.Pt();
313
314           if (c->Chi2() == 100) // MC particle
315             mcpt += nP.Pt();
316
317           jet->AddClusterAt(cid, nc);
318           ++nc;
319           ++nneutral;
320         }
321       } else {
322         AliError(Form("%s: No logical way to end up here.", GetName()));
323         continue;
324       }
325     }
326     jet->SetNumberOfTracks(nt);
327     jet->SetNumberOfClusters(nc);
328     jet->SortConstituents();
329     jet->SetMaxNeutralPt(maxNe);
330     jet->SetMaxChargedPt(maxCh);
331     jet->SetNEF(neutralE / jet->E());
332     jet->SetArea(fjw.GetJetArea(ij));
333     jet->SetNumberOfCharged(ncharged);
334     jet->SetNumberOfNeutrals(nneutral);
335     jet->SetMCPt(mcpt);
336     if (gall > 0)
337       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
338     else 
339       jet->SetAreaEmc(-1);
340     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
341         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
342         (jet->Eta() > geom->GetArm1EtaMin()) && 
343         (jet->Eta() < geom->GetArm1EtaMax()))
344       jet->SetAxisInEmcal(kTRUE);
345     jetCount++;
346   }
347 }