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