setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskEmcalJetSample.cxx
CommitLineData
8628b70c 1// $Id$
9f52c61f 2//
3// Jet sample analysis task.
4//
5// Author: S.Aiola
6
7#include <TClonesArray.h>
8#include <TH1F.h>
9#include <TH2F.h>
0925db21 10#include <TH3F.h>
9f52c61f 11#include <TList.h>
12#include <TLorentzVector.h>
13
14#include "AliVCluster.h"
4bb8c6df 15#include "AliAODCaloCluster.h"
0925db21 16#include "AliESDCaloCluster.h"
9f52c61f 17#include "AliVTrack.h"
18#include "AliEmcalJet.h"
19#include "AliRhoParameter.h"
20#include "AliLog.h"
9239b066 21#include "AliJetContainer.h"
8612dfc8 22#include "AliParticleContainer.h"
23#include "AliClusterContainer.h"
0925db21 24#include "AliPicoTrack.h"
9f52c61f 25
26#include "AliAnalysisTaskEmcalJetSample.h"
27
28ClassImp(AliAnalysisTaskEmcalJetSample)
29
30//________________________________________________________________________
31AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
8612dfc8 32 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
4bb8c6df 33 fHistTracksPt(0),
34 fHistClustersPt(0),
35 fHistLeadingJetPt(0),
36 fHistJetsPhiEta(0),
37 fHistJetsPtArea(0),
38 fHistJetsPtLeadHad(0),
39 fHistJetsCorrPtArea(0),
0925db21 40 fHistPtDEtaDPhiTrackClus(0),
41 fHistPtDEtaDPhiClusTrack(0),
8612dfc8 42 fJetsCont(0),
43 fTracksCont(0),
44 fCaloClustersCont(0)
9f52c61f 45
46{
47 // Default constructor.
48
4bb8c6df 49 fHistTracksPt = new TH1*[fNcentBins];
50 fHistClustersPt = new TH1*[fNcentBins];
51 fHistLeadingJetPt = new TH1*[fNcentBins];
52 fHistJetsPhiEta = new TH2*[fNcentBins];
53 fHistJetsPtArea = new TH2*[fNcentBins];
54 fHistJetsPtLeadHad = new TH2*[fNcentBins];
55 fHistJetsCorrPtArea = new TH2*[fNcentBins];
56
57 for (Int_t i = 0; i < fNcentBins; i++) {
9f52c61f 58 fHistTracksPt[i] = 0;
59 fHistClustersPt[i] = 0;
60 fHistLeadingJetPt[i] = 0;
61 fHistJetsPhiEta[i] = 0;
62 fHistJetsPtArea[i] = 0;
63 fHistJetsPtLeadHad[i] = 0;
64 fHistJetsCorrPtArea[i] = 0;
65 }
66
67 SetMakeGeneralHistograms(kTRUE);
68}
69
70//________________________________________________________________________
71AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
8612dfc8 72 AliAnalysisTaskEmcalJet(name, kTRUE),
4bb8c6df 73 fHistTracksPt(0),
74 fHistClustersPt(0),
75 fHistLeadingJetPt(0),
76 fHistJetsPhiEta(0),
77 fHistJetsPtArea(0),
78 fHistJetsPtLeadHad(0),
79 fHistJetsCorrPtArea(0),
0925db21 80 fHistPtDEtaDPhiTrackClus(0),
81 fHistPtDEtaDPhiClusTrack(0),
8612dfc8 82 fJetsCont(0),
83 fTracksCont(0),
84 fCaloClustersCont(0)
9f52c61f 85{
86 // Standard constructor.
87
4bb8c6df 88 fHistTracksPt = new TH1*[fNcentBins];
89 fHistClustersPt = new TH1*[fNcentBins];
90 fHistLeadingJetPt = new TH1*[fNcentBins];
91 fHistJetsPhiEta = new TH2*[fNcentBins];
92 fHistJetsPtArea = new TH2*[fNcentBins];
93 fHistJetsPtLeadHad = new TH2*[fNcentBins];
94 fHistJetsCorrPtArea = new TH2*[fNcentBins];
9239b066 95
4bb8c6df 96 for (Int_t i = 0; i < fNcentBins; i++) {
9f52c61f 97 fHistTracksPt[i] = 0;
98 fHistClustersPt[i] = 0;
99 fHistLeadingJetPt[i] = 0;
100 fHistJetsPhiEta[i] = 0;
101 fHistJetsPtArea[i] = 0;
102 fHistJetsPtLeadHad[i] = 0;
103 fHistJetsCorrPtArea[i] = 0;
104 }
105
106 SetMakeGeneralHistograms(kTRUE);
107}
108
109//________________________________________________________________________
110AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
111{
112 // Destructor.
113}
114
115//________________________________________________________________________
116void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
117{
118 // Create user output.
119
120 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
121
4bb8c6df 122 fJetsCont = GetJetContainer(0);
123 if(fJetsCont) { //get particles and clusters connected to jets
124 fTracksCont = fJetsCont->GetParticleContainer();
125 fCaloClustersCont = fJetsCont->GetClusterContainer();
126 } else { //no jets, just analysis tracks and clusters
127 fTracksCont = GetParticleContainer(0);
128 fCaloClustersCont = GetClusterContainer(0);
129 }
130 fTracksCont->SetClassName("AliVTrack");
05f22af9 131 fCaloClustersCont->SetClassName("AliVCluster");
8612dfc8 132
9f52c61f 133 TString histname;
134
4bb8c6df 135 for (Int_t i = 0; i < fNcentBins; i++) {
9239b066 136 if (fParticleCollArray.GetEntriesFast()>0) {
9f52c61f 137 histname = "fHistTracksPt_";
138 histname += i;
139 fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
140 fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
141 fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
142 fOutput->Add(fHistTracksPt[i]);
143 }
144
9239b066 145 if (fClusterCollArray.GetEntriesFast()>0) {
9f52c61f 146 histname = "fHistClustersPt_";
147 histname += i;
148 fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
149 fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
150 fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
151 fOutput->Add(fHistClustersPt[i]);
152 }
153
9239b066 154 if (fJetCollArray.GetEntriesFast()>0) {
9f52c61f 155 histname = "fHistLeadingJetPt_";
156 histname += i;
157 fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
158 fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
159 fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
160 fOutput->Add(fHistLeadingJetPt[i]);
161
162 histname = "fHistJetsPhiEta_";
163 histname += i;
164 fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
165 fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
166 fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
167 fOutput->Add(fHistJetsPhiEta[i]);
168
169 histname = "fHistJetsPtArea_";
170 histname += i;
171 fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
172 fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
173 fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
174 fOutput->Add(fHistJetsPtArea[i]);
175
176 histname = "fHistJetsPtLeadHad_";
177 histname += i;
178 fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
179 fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
180 fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
181 fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
182 fOutput->Add(fHistJetsPtLeadHad[i]);
183
9239b066 184 if (!(GetJetContainer()->GetRhoName().IsNull())) {
9f52c61f 185 histname = "fHistJetsCorrPtArea_";
186 histname += i;
187 fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
188 fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
189 fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
190 fOutput->Add(fHistJetsCorrPtArea[i]);
191 }
192 }
193 }
0925db21 194
195 histname = "fHistPtDEtaDPhiTrackClus";
196 fHistPtDEtaDPhiTrackClus = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
197 fOutput->Add(fHistPtDEtaDPhiTrackClus);
198
199 histname = "fHistPtDEtaDPhiClusTrack";
200 fHistPtDEtaDPhiClusTrack = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
201 fOutput->Add(fHistPtDEtaDPhiClusTrack);
202
9f52c61f 203 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
204}
205
206//________________________________________________________________________
207Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
208{
209 // Fill histograms.
210
8612dfc8 211 if (fTracksCont) {
4bb8c6df 212 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
8612dfc8 213 while(track) {
9f52c61f 214 fHistTracksPt[fCentBin]->Fill(track->Pt());
4bb8c6df 215 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
9f52c61f 216 }
217 }
218
8612dfc8 219 if (fCaloClustersCont) {
0925db21 220 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
8612dfc8 221 while(cluster) {
9f52c61f 222 TLorentzVector nPart;
223 cluster->GetMomentum(nPart, fVertex);
224 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
0925db21 225
226 cluster = fCaloClustersCont->GetNextAcceptCluster();
9f52c61f 227 }
228 }
229
8612dfc8 230 if (fJetsCont) {
231 AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
232 while(jet) {
9239b066 233
9f52c61f 234 fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
235 fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
236
22339f55 237 Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
9f52c61f 238 fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
239
388f9cd5 240 if (fHistJetsCorrPtArea[fCentBin]) {
241 Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
242 fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
243 }
8612dfc8 244 jet = fJetsCont->GetNextAcceptJet();
9f52c61f 245 }
8612dfc8 246
247 jet = fJetsCont->GetLeadingJet();
6b840003 248 if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
9f52c61f 249 }
250
0925db21 251 CheckClusTrackMatching();
252
9f52c61f 253 return kTRUE;
254}
255
256//________________________________________________________________________
0925db21 257void AliAnalysisTaskEmcalJetSample::CheckClusTrackMatching()
258{
259
260 if(!fTracksCont || !fCaloClustersCont)
261 return;
262
263 Double_t deta = 999;
264 Double_t dphi = 999;
265
266 //Get closest cluster to track
267 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
268 while(track) {
269 //Get matched cluster
270 Int_t emc1 = track->GetEMCALcluster();
271 if(fCaloClustersCont && emc1>=0) {
272 AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
273 if(clusMatch) {
274 AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
275 fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
276 }
277 }
278 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
279 }
280
281 //Get closest track to cluster
282 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
283 while(cluster) {
284 TLorentzVector nPart;
285 cluster->GetMomentum(nPart, fVertex);
286 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
287
288 //Get matched track
289 AliVTrack *mt = NULL;
290 AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
291 if(acl) {
292 if(acl->GetNTracksMatched()>1)
293 mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
294 }
295 else {
296 AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
297 Int_t im = ecl->GetTrackMatchedIndex();
298 if(fTracksCont && im>=0) {
299 mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
300 }
301 }
302 if(mt) {
303 AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
304 fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
305
306 /* //debugging
307 if(mt->IsEMCAL()) {
308 Int_t emc1 = mt->GetEMCALcluster();
309 Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
310 AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
311 AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
312 Printf("deta: %f dphi: %f",deta,dphi);
313 }
314 */
315 }
316 cluster = fCaloClustersCont->GetNextAcceptCluster();
317 }
318}
319
320//________________________________________________________________________
8612dfc8 321void AliAnalysisTaskEmcalJetSample::ExecOnce() {
322
323 AliAnalysisTaskEmcalJet::ExecOnce();
324
325 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
326 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
327 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
328
329}
330
331//________________________________________________________________________
9f52c61f 332Bool_t AliAnalysisTaskEmcalJetSample::Run()
333{
334 // Run analysis code here, if needed. It will be executed before FillHistograms().
335
336 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
337}
338
339//________________________________________________________________________
340void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *)
341{
342 // Called once at the end of the analysis.
343}