Cosmetics
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliEmcalJetTask.cxx
CommitLineData
48d874ff 1// $Id$
542ac102 2//
3// Emcal jet finder task.
4//
5//
48d874ff 6
914d486c 7#include "AliEmcalJetTask.h"
6ea168d0 8#include <TChain.h>
48d874ff 9#include <TClonesArray.h>
6ea168d0 10#include <TH1F.h>
11#include <TH2F.h>
12#include <TList.h>
13#include <TLorentzVector.h>
48d874ff 14#include <TParticle.h>
48d874ff 15#include "AliAnalysisManager.h"
6ea168d0 16#include "AliCentrality.h"
35789a2d 17#include "AliVCluster.h"
18#include "AliVTrack.h"
914d486c 19#include "AliEmcalJet.h"
48d874ff 20#include "AliFJWrapper.h"
48d874ff 21
914d486c 22ClassImp(AliEmcalJetTask)
48d874ff 23
48d874ff 24//________________________________________________________________________
914d486c 25AliEmcalJetTask::AliEmcalJetTask() :
26 AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0 27 fTracksName("Tracks"),
28 fCaloName("CaloClusters"),
29 fJetsName("Jets"),
30 fAlgo(1),
31 fRadius(0.4),
32 fType(0),
914d486c 33 fMinJetTrackPt(0.15),
34 fMinJetClusPt(0.15),
6ea168d0 35 fJets(0)
36{
37 // Default constructor.
38
39 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
40}
41
42//________________________________________________________________________
914d486c 43AliEmcalJetTask::AliEmcalJetTask(const char *name) :
44 AliAnalysisTaskSE("AliEmcalJetTask"),
7efbea04 45 fTracksName("Tracks"),
48d874ff 46 fCaloName("CaloClusters"),
7efbea04 47 fJetsName("Jets"),
48d874ff 48 fAlgo(1),
49 fRadius(0.4),
50 fType(0),
6ea168d0 51 fMinJetTrackPt(0.15),
914d486c 52 fMinJetClusPt(0.15),
7efbea04 53 fJets(0)
48d874ff 54{
55 // Standard constructor.
7efbea04 56
48d874ff 57 if (!name)
58 return;
7efbea04 59
48d874ff 60 SetName(name);
7efbea04 61 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
6ea168d0 62
27fde6c7 63 //DefineInput(0,TChain::Class());
64 //DefineOutput(1,TList::Class());
48d874ff 65}
66
67//________________________________________________________________________
914d486c 68AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 69{
70 // Destructor
71}
72
73//________________________________________________________________________
914d486c 74void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff 75{
76 // Create user objects.
77
914d486c 78 fJets = new TClonesArray("AliEmcalJet");
48d874ff 79 fJets->SetName(fJetsName);
48d874ff 80}
81
82//________________________________________________________________________
914d486c 83void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 84{
85 // Main loop, called for each event.
7efbea04 86 // Add jets to event if not yet there
48d874ff 87
48d874ff 88 if (!(InputEvent()->FindListObject(fJetsName)))
89 InputEvent()->AddObject(fJets);
51914326 90 else {
91 // IMPORTANT: if it is not an AliESDEvent, non-std TClonesArray will not be cleared automatically!
92 if (!(InputEvent()->InheritsFrom("AliESDEvent")))
93 fJets->Delete();
94 }
48d874ff 95
96 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
97 TClonesArray *tracks = 0;
98 TClonesArray *clus = 0;
99 TList *l = InputEvent()->GetList();
6ea168d0 100
101 Float_t cent=100;
102 AliCentrality *centrality = dynamic_cast<AliCentrality*>(l->FindObject("Centrality"));
103 if(centrality)
104 cent = centrality->GetCentralityPercentile("V0M");
105
48d874ff 106 if ((fType==0)||(fType==1)) {
7efbea04 107 if (fTracksName == "Tracks")
48d874ff 108 am->LoadBranch("Tracks");
7efbea04 109 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
48d874ff 110 if (!tracks) {
7efbea04 111 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
48d874ff 112 return;
113 }
114 }
115 if ((fType==0)||(fType==2)) {
116 if (fCaloName == "CaloClusters")
117 am->LoadBranch("CaloClusters");
118 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
119 if (!clus) {
120 AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
121 return;
122 }
48d874ff 123 }
124
6ea168d0 125 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 126}
127
128//________________________________________________________________________
914d486c 129void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 130{
131 // Called once at the end of the analysis.
48d874ff 132}
133
134//________________________________________________________________________
914d486c 135void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 136{
137 // Find jets.
138
139 TString name("kt");
140 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
141 if (algo>=1) {
142 name = "antikt";
143 jalgo = fastjet::antikt_algorithm;
144 }
145
146 AliFJWrapper fjw(name, name);
147 fjw.SetR(radius);
6ea168d0 148 fjw.SetAlgorithm(jalgo);
48d874ff 149 fjw.SetMaxRap(0.9);
150 fjw.Clear();
151
152 if (tracks) {
153 const Int_t Ntracks = tracks->GetEntries();
154 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
7efbea04 155 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
48d874ff 156 if (!t)
157 continue;
914d486c 158 if (t->Pt()<fMinJetTrackPt)
159 continue;
48d874ff 160 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
161 }
162 }
6ea168d0 163
48d874ff 164 if (clus) {
165 Double_t vertex[3] = {0, 0, 0};
166 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
167 const Int_t Nclus = clus->GetEntries();
914d486c 168 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 169 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
48d874ff 170 if (!c->IsEMCAL())
171 continue;
172 TLorentzVector nPart;
173 c->GetMomentum(nPart, vertex);
174 Double_t energy = nPart.P();
914d486c 175 if (energy<fMinJetClusPt)
176 continue;
6ea168d0 177 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), energy, -iClus-1);
48d874ff 178 }
179 }
180
7efbea04 181 // run jet finder
48d874ff 182 fjw.Run();
6ea168d0 183
48d874ff 184 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
185 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
186 if (jets_incl[ij].perp()<1)
187 continue;
914d486c 188 AliEmcalJet *jet = new ((*fJets)[jetCount])
189 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
900f5135 190 jet->SetArea(fjw.GetJetArea(ij));
6ea168d0 191 Double_t vertex[3] = {0, 0, 0};
192 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 193 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
6ea168d0 194 Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
35789a2d 195 jet->SetNumberOfTracks(constituents.size());
196 jet->SetNumberOfClusters(constituents.size());
197 Int_t nt = 0;
198 Int_t nc = 0;
48d874ff 199 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
200 Int_t uid = constituents[ic].user_index();
6ea168d0 201 if (uid>=0){
35789a2d 202 jet->AddTrackAt(uid, nt);
6ea168d0 203 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
204 if (t->Pt()>maxTrack)
205 maxTrack=t->Pt();
35789a2d 206 nt++;
6ea168d0 207 } else {
35789a2d 208 jet->AddClusterAt(-(uid+1),nc);
6ea168d0 209 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(-(uid+1)));
210 TLorentzVector nP;
211 c->GetMomentum(nP, vertex);
212 neutralE += nP.P();
213 if (nP.P()>maxCluster)
214 maxCluster=nP.P();
35789a2d 215 nc++;
48d874ff 216 }
217 }
35789a2d 218 jet->SetNumberOfTracks(nt);
219 jet->SetNumberOfClusters(nc);
6ea168d0 220 jet->SetMaxTrackPt(maxTrack);
221 jet->SetMaxClusterPt(maxCluster);
48d874ff 222 jet->SetNEF(neutralE/jet->E());
48d874ff 223 jetCount++;
224 }
225}