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