Changes for #93556: Trigger list in Pb-Pb 2011 simulation
[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"
13#include "AliESDCaloCluster.h"
48d874ff 14#include "AliESDtrack.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);
86
87 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
88 TClonesArray *tracks = 0;
89 TClonesArray *clus = 0;
90 TList *l = InputEvent()->GetList();
6ea168d0 91
92 Float_t cent=100;
93 AliCentrality *centrality = dynamic_cast<AliCentrality*>(l->FindObject("Centrality"));
94 if(centrality)
95 cent = centrality->GetCentralityPercentile("V0M");
96
48d874ff 97 if ((fType==0)||(fType==1)) {
7efbea04 98 if (fTracksName == "Tracks")
48d874ff 99 am->LoadBranch("Tracks");
7efbea04 100 tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
48d874ff 101 if (!tracks) {
7efbea04 102 AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
48d874ff 103 return;
104 }
105 }
106 if ((fType==0)||(fType==2)) {
107 if (fCaloName == "CaloClusters")
108 am->LoadBranch("CaloClusters");
109 clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
110 if (!clus) {
111 AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
112 return;
113 }
48d874ff 114 }
115
6ea168d0 116 FindJets(tracks, clus, fAlgo, fRadius, cent);
48d874ff 117}
118
119//________________________________________________________________________
914d486c 120void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 121{
122 // Called once at the end of the analysis.
48d874ff 123}
124
125//________________________________________________________________________
914d486c 126void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
48d874ff 127{
128 // Find jets.
129
130 TString name("kt");
131 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
132 if (algo>=1) {
133 name = "antikt";
134 jalgo = fastjet::antikt_algorithm;
135 }
136
137 AliFJWrapper fjw(name, name);
138 fjw.SetR(radius);
6ea168d0 139 fjw.SetAlgorithm(jalgo);
48d874ff 140 fjw.SetMaxRap(0.9);
141 fjw.Clear();
142
143 if (tracks) {
144 const Int_t Ntracks = tracks->GetEntries();
145 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
7efbea04 146 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
48d874ff 147 if (!t)
148 continue;
914d486c 149 if (t->Pt()<fMinJetTrackPt)
150 continue;
48d874ff 151 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
152 }
153 }
6ea168d0 154
48d874ff 155 if (clus) {
156 Double_t vertex[3] = {0, 0, 0};
157 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
158 const Int_t Nclus = clus->GetEntries();
914d486c 159 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
7efbea04 160 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
48d874ff 161 if (!c->IsEMCAL())
162 continue;
163 TLorentzVector nPart;
164 c->GetMomentum(nPart, vertex);
165 Double_t energy = nPart.P();
914d486c 166 if (energy<fMinJetClusPt)
167 continue;
6ea168d0 168 fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), energy, -iClus-1);
48d874ff 169 }
170 }
171
7efbea04 172 // run jet finder
48d874ff 173 fjw.Run();
6ea168d0 174
48d874ff 175 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
176 for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
177 if (jets_incl[ij].perp()<1)
178 continue;
914d486c 179 AliEmcalJet *jet = new ((*fJets)[jetCount])
180 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
900f5135 181 jet->SetArea(fjw.GetJetArea(ij));
6ea168d0 182 Double_t vertex[3] = {0, 0, 0};
183 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
48d874ff 184 vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
6ea168d0 185 Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
48d874ff 186 for(UInt_t ic=0; ic<constituents.size(); ++ic) {
187 Int_t uid = constituents[ic].user_index();
6ea168d0 188 if (uid>=0){
189 AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
190 if (t->Pt()>maxTrack)
191 maxTrack=t->Pt();
192 } else {
193 AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(-(uid+1)));
194 TLorentzVector nP;
195 c->GetMomentum(nP, vertex);
196 neutralE += nP.P();
197 if (nP.P()>maxCluster)
198 maxCluster=nP.P();
48d874ff 199 }
200 }
6ea168d0 201 jet->SetMaxTrackPt(maxTrack);
202 jet->SetMaxClusterPt(maxCluster);
48d874ff 203 jet->SetNEF(neutralE/jet->E());
48d874ff 204 jetCount++;
205 }
206}